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The stability of a stratified flow 


By WILLIAM PAUL GRAEBEL 


Department of Engineering Mechanics, University of Michigan, Ann Arbor, Michigan 
(Received 20 November 1959) 


This paper deals theoretically with the problem of the hydrodynamic stability 
of a stratified flow of a viscous fluid. The primary flow consists of two laminar 
streams of viscous fluids of different densities flowing in opposite directions 
between two parallel inclined planes under the action of gravity. The effect of 
surface tension at the interface of the two fluids is included in the formulation 
of the problem. 

Since instability can be expected to occur at low Reynolds numbers when the 
inclination is nearly vertical, the solution of the Orr-Sommerfeld equations is 
developed as a power series in the transverse space co-ordinate. It is shown that 
for the vertical case, the flow is unstable for all values of the Reynolds number. 
Surface tension is found to influence both the direction and celerity of the dis- 
turbance. Results are also given for inclinations slightly away from the vertical, 
where small critical Reynolds numbers do exist. 


1. Introduction 

For most parallel flows past fixed boundaries whose stability has been in- 
vestigated so far, instability occurs at rather large Reynolds numbers. The Orr— 
Sommerfeld equation (see Lin 1955) governing the stability of such flows must 
therefore be solved for large values of the Reynolds number, which appears as 
a parameter in these equations. The asymptotic solutions of the Orr-Sommerfeld 
equation, appropriate for large Reynolds numbers, have singularities at the 
critical points where the wave velocity of the disturbance is equal to the velocity 
of the mean flow. Great care must therefore be exercised in the evaluation of 
these solutions as a critical point is crossed. These singularities of the solutions 
are, however, not inherent in the Orr-Sommerfeld equation, and are introduced 
entirely by the method of solution. For flows which can be expected to be 
unstable at low Reynolds numbers, the appropriate solutions can be expressed 
in ascending powers of the Reynolds number or of one of the co-ordinates. Since 
asymptotic solutions are not needed, the aforementioned singularities do not 
occur. The study of the stability of flows which can be expected to become 
unstable at low Reynolds numbers can therefore be carried out by conventional 
methods. In view of this, it is perhaps somewhat surprising that until recent 
years problems of hydrodynamic instability at low Reynolds numbers have been 
neglected by research workers. 

The flow whose stability is studied here is a stratified flow of two fluids of 


equal viscosity but different densities. It is entirely motivated by gravity and 
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its velocity distribution is antisymmetric, with the lighter fluid flowing in a 
direction opposite to that of the heavier fluid, with a point of inflexion of the 
velocity profile at the interface. Since a point of inflexion is known to have a 
destabilizing effect, and since the work of Yih (1954) and of Benjamin (1957) 
has shown that a free surface also has a destabilizing effect, it can be expected 
to be unstable at low Reynolds numbers, especially if the slope of the mean flow 
is steep. The Orr-Sommerfeld equation will therefore be solved for small Reynolds 
numbers. 

The type of flow investigated here is encountered in extraction columns in 
chemical engineering, and can also be found in closed channels such as tunnels 
and mine shafts. It can also be expected that the results have some bearing on 
stratified flows occurring in the oceans and the atmosphere, where layers of 
warm water or air flow over and counter to colder water or air (i.e. the flow of air 
after cold fronts). 


2. The primary flow 


The present investigation concerns the stability of a steady laminar stratified 
flow of an incompressible viscous fluid between two parallel fixed planes. The 
spacing of the planes is denoted by 2b. The origin of the co-ordinate system is 
taken half-way between the planes, with the X-axis parallel to the planes. The 
planes are inclined at an angle @ with the horizontal (see figure 1). 





FIGURE 1. Diagram of the primary flow, showing the co-ordinate 
axes and the velocity profile. 


The fluid occupying the region 0 < Y < b, with density p,, flows up the inclined 
plane in the direction of negative X. The fluid occupying the region —b < Y < 0, 
with density p, greater than p,, flows down the plane under the action of gravity. 
The viscosity ~ of the two fluids are considered equal. Gravity is the sole 
motivating force for the flow, with the heavier fluid displacing the lighter fluid 
in a reservoir at X = +00. The volumetric discharge across the channel is 
therefore zero. 
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The stability of a stratified flow 


The Navier-Stokes equations which govern the primary flow are 


0=—-P, x+p,9sind+W, yy, ) 

0=-—P,;-p,gcosO for 0< Y<b,J ) 
and 0=—-—P, y+pogsind+pnU, yy, (2) 

0=—-P,y—-p.gceos@ for —-b< Y 


where U, and U,, the components of velocity of the two fluids in the X-direction, 
are functions of Y only, and letter subscripts following commas denote partial 
differentiation with respect to that quantity. The components of velocity in the 
Y- and Z-directions are zero. The condition of continuity is thus automatically 
satisfied. Pressure is indicated by P and the gravitational acceleration by g. 
These equations are easily solved for U, and U,. It can easily be verified that 





0, =U, = ; <j) - éren, (3) 

1 U, = —U. : b<Y<0 4 
: ,- wf T+Q]] crereo “ 
where Uy = (Pe af ) 52 2g sin 0 (5) 


j 


is four times the maximum velocity. The boundary conditions corresponding 
to (3) and (4) are conditions of no slipping at the fixed boundaries. The condition 
of zero volumetric discharge is also satisfied by (3) and (4). The pressure gradient 
in the X-direction is 

Pix = Phx = HMp2tpy)gsind. 


Thus, the X-variation of pressure is the same as that of the hydrostatic pressure 
in a fluid of density 4(p,+ ,), while the Y-variation corresponds to a hydrostatic 
pressure in the fluids at present under consideration. 


3. Formulation of the stability problem 
A. Differential equations 

The basic equations to be satisfied are the Navier-Stokes equations and the 
continuity equations. Following the usual approach to stability problems, it is 
assumed that the perturbed velocity and pressure fields are expressible in power 
series expansions in terms of an amplitude parameter ¢, assumed small and 
constant (see Lin 1955, § 1.2). Upon substituting these expansions into the field 
equations and setting coefficients of powers of ¢€ to zero, it is noted that the 
equations in ¢° are those already described as governing the primary flow. The 
equations in e! are 


UF 7 + U, Wx + U, yu* Pt x ee Vays ; (6) 


y+ LVet 
Pi 
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where the disturbance is assumed to be a two-dimensional one. The convention 
has been adopted that capital letters refer to the primary flow and small letters 
with an asterisk in the upper right-hand corner refer to the first-order disturbance 
flow. The equations hold in the upper region when 7 = 1 and in the lower region 
when i = 2. As is customary in stability analyses, attention will be confined 
solely to the first-order equations, with the assumption that for sufficiently small 
disturbances the linearized equations can be used without appreciable error. 
To facilitate the solution of (6), (7) and (8), it is convenient to put them in 
dimensionless form. The following dimensionless parameters are introduced 


us ve pi m X ; r 
ue=—, V,=—, D,= ee = —_, 
‘  Ony ' Uy ™ PiU F's | ( 
pea ad ati 9) 
ee, Gest, Bae, 4o8. | 
b Uy lt Pi 


Formulating (6), (7) and (8) in terms of these parameters, the following dimen- 
sionless equations are obtained 


: : ] 

UW; pt l iVing l iyi i Piz ~ R V2u;, (10) 
7m 

Vit + C Ving seat Piy ? RY (11) 

U; t+ Viy = @, (12) 


in which R; is the Reynolds number based on the half-spacing 6 and Uj,. Since 
the flow is motivated by gravity, it is expected that the Reynolds number and 
the Froude number F are not independent of one another. Such is indeed the 
case, for if F is defined by F? = U4,/gb, then from the primary flow (5) provides 
the relationship F? = i(r—1)R,sin0. (13) 

To reduce the number of equations which are to be solved, it is convenient to 
introduce stream functions. To simplify solution of (10), (11) and (12), the 
solutions will be assumed to have exponential time factors. Further, to reduce 
the partial differential equations to ordinary differential equations, it is assumed 
that the disturbance may be resolved into Fourier components, which are thus 
periodic in x. A general disturbance then would consist of either a Fourier series 
or a Fourier integral of such components. To determine the stability of the flow 
it is sufficient to consider the effect of a single Fourier component of general 
period. 

Only two-dimensional disturbances are considered, since the work of Squire 
(1933), Yih (1955) and others has shown that the stability or instability of a 
three-dimensional disturbance can be determined from that of a two-dimen- 
sional disturbance at a higher Reynolds number. 

The velocity components and the pressure can thus be assumed to have the 
following forms, in which f(y) and h(y) are introduced as ‘stream functions’ to 
satisfy (12) automatically ; 

. “ Uy = f’(y) etez—c, | 
ty = ~iafly)eiae-,| 


q(y) etalz—et) 


(14) 
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P1 





























The stability of a stratified flow 
in0O<y<1,and 
Uy = h'(y) eta(e—el) | 


V2 = —tah(y) eta, & (15) 


Pz = Qoly) er 


in —1 <y <0, where primes denote differentiation with respect to y. The 
dimensionless wave number « is related to the wave length A of the disturbance 
by « = 27b/A. The real part of c represents the dimensionless disturbance celerity, 
and a times the imaginary part of c gives the dimensionless growth rate of the 
disturbance. 

Upon solving (10), (11), (14) and (15) for the pressure terms, the following 
are obtained 1 , 
= gp S"— at) +(e— Ups’ + Uf, 

ak, 


(16) 
Re ar” _ yar 2(n _ TT 
1 iar, 4 atf)+a%(c—U,)f | 
in0O<y<1,and 
] 
q2 = ial, (h” —a*h’) + (e—U,)h'+ ad sa 
7 


’ l oL" 9 ry 
q2 = ma —ath)+a?(e—U,)h | 


in —1 < y < 0. The pressure terms can now be eliminated to obtain 
fi¥+[—2a%+iaR,(c—U,)]f" + [at +iaR, Ul +108 R,(U,—c)] f = 0 (18) 
inO<y<1,and 
hiv + [—20?+iaR,(c—U,)]h" + [at+iaR, UZ+iaeR(U,—c)]h=9 (19) 
in —1 <y <0. These equations are the Orr-Sommerfeld equations. From 
(3), (4) and (9) the functions U, and U, are 


U,=—-yt+y, U,=-y-y’. (20) 


B. Boundary conditions 


Since the differential equations (18) and (19) are two in number and each is 
of the fourth order, there must be eight boundary conditions imposed to specify 
the mathematical problem completely. Two conditions are imposed at each of 
the fixed boundaries, and a total of four are imposed at the interface. 

Since the fluids are considered viscous, there must be no slip at the fixed 
boundaries. Hence, in terms of the stream functions, 


fy =9, f/U)=9, A(-1)=0, h'(-1)=0. (21) 


The interface is assumed to be displaced from the z-axis by a small amount 
which in dimensional form is denoted by by. This introduces a further unknown 
into the problem, which, however, may be readily determined from the kine- 
matical condition that the component of velocity at the interface must equal the 
time derivative of 7. In what follows, quantities evaluated at the interface are 
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expanded in a Taylors series about y = 0, and 7 is considered to be of the same 
order as the velocity disturbances. In keeping with the previous linearization, 
only first-order disturbance terms are retained. 

At the interface, the velocity components must be continuous; hence 


(0) = A(9), f’(0) = A'(0). (22) 
The kinematic condition at the interface is, in dimensionless terms, 
deel (23) 
at y = 0. Then N= —vaf(0) ete 
and upon solving for 7 
n = _ fo) ei, (24) 


The shear stress must also be continuous across the interface. To the first 
order, the dimensionless shear stress at the interface is given by 


Try = a(U" +u'+v ,+U"y). (25) 
t 
Therefore the boundary condition imposed is, after some simplification, 
4f(0) + ¢f"(0)—ch"(0) = 0. (26) 
The dimensionless normal stress at the interface is 


l 


z. — — 
pus, 


> pod a. Or 
yy — (P+P yby + 2uU yr") — p, (27) 
where P is the pressure in the primary flow. The condition imposed is that the 
difference in normal stresses must be equal to the curvature of the interface times 
y, the surface tension. Since P must already be continuous, and since 


Py =—gpcos#, 


the boundary condition becomes 


bgcosO a? 


l[(r-1 
91(9) — r92(0) + - i +p f(0) = 9, (28) 


72 
M 


where W is the Weber number defined by W = p,bU%,/y. This can be written 
entirely in terms of the stream functions by using (16) and (17). The resulting 
equation is, after further simplification, 


c[f"(0) —h"(0)] + (r— 1) ia Ry e[f(0) — ef (0)] + ia[4 cot 6 + 22S] f(0) = 0. (29) 


The parameter S is defined by 


R, y 
y= = —,, 30 
W wy (30) 


and represents the ratio of surface tension to viscous forces. 














fee - 
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4. Solution of the stability problem 

A. General 
For the purpose of solving (18) and (19), much labour can be saved by employing 
the following notation: 





; ( f(y) in O<y<l, 
Ply) = | ; 
lLa(y) in —l<y<0, 
B= a? 
4 (“ =iaR, in O<y<l, ;¢ (31) 
~  (4,=iaR, in -l<y<0, 
POM i n O<y<l 
i 1 in —l<y<0 


Then (18) and (19) can both be solved simultaneously by solving 
gi¥+[—-2B+ A(c+y—dy?)] 6" + [B2+2Ad6+ AB(—c—y+ dy*)]¢ = 90. (32) 


As discussed in the Introduction, instability can be expected to occur at low 
Reynolds numbers when the inclination 6 approaches 90°. Two techniques at 
once are suggested as alternate methods of solution. The first, expansion of ¢ 
in ascending powers of the Reynolds number, is certainly the more typical of 
stability problems and was used by Yih (1954). The alternate approach, expan- 
sion of ¢ in powers of the co-ordinate y, has, to the best of the author’s know- 
ledge, been previously employed successfully only by Benjamin (1957) in spite 
of its early recognition by Kelvin (1887) and others. This method seems more 
useful in carrying out the solution for ¢ to higher powers of the Reynolds number, 
and conclusions are drawn on the basis of these calculations. Of course, the 
validity of approximations based on either method can be determined only 
a posteriori. 

B. Expansion in powers of the co-ordinate y 


Following Benjamin (1957), a solution of the form 
oy) = X any" (33) 


can be assumed. It is most convenient to expand around the origin since four 
of the boundary conditions must be satisfied there. Upon substituting (33) into 
(32) and setting coefficients of y” to zero, the following recurrence relation is 
obtained: 

(n—4)! : ees ; 
— {(n — 2) (n—3) (2B—Ac)a,_.—(n—4) (n—3) Aa, _, 


n ' 


nN: 


+ [Ad(n —4) (n—5)— B?-2Ad+ABc]a,_,+ABa,_,—ABéa,_,}. (34) 


All of the a,, for n > 4 can now be found from (34) in terms of dp, a, a, and a3. 
Since these four coefficients are arbitrary, four independent solutions of the 
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fourth-order equation are thus obtained. They are (where do, @,, 4, and a, have 
been set equal to 1, 1, 1 and 1/3!, respectively) 
,(y) = 1+ (1/4!) (A Be — B?— 2Ad) y* + (1/5!) ABy? 
+ (1/6!) (2A%ed —- 6A BS + 3A Bee — A* Be — 2B) 76 
+ (1/7!) (6426+ 5A B? —4A?Be?) y? 
+ (1/8!) (—20A?—3B!+ 64 B8c — 44? B2c? — 20A B76 + 18A?2Bed 
— 2A%c?$ + A Bc3 — 4A?B) y8 
+(1/9!) (144 B3 + 5042B6 + 9A3Be? — 22.42B2c — 1643cd) y9 +..., (35) 
d(y) = y+ (1/5!) (ABe— B2—2A0) y> + (1/6!) 2ABy® 
+ (1/7!) (2A?ed — 10A Bd + 34 Bee — A? Bc? — 2 B3) y? 
+ (1/8!) (8426+ 8A B?-6A?Be) 
+ (1/9!) (—36A?— 3B! + 6A Bc — 364A B26 + 830A2Bed — 4A? B?c? 
— 2A%c?6 — 10A?B + A Bc?) y?® 
+ (1/10!) (204 B8 + 122A?BS + 1242 Be? — 30A2B2e — 20A3cd) y! + ..., 
(36) 
3(y) = y? + (1/4!) (4B — 2Ac) y*— (1/5!) 2A y® 
+ (1/6!) (6B? —6A Be + 2A?c?) y® + (1/7!) (8A%ce — LOA B) y’ 
+ (1/8!) (8A? + 164 Bd — 20A2cd + 8B3 — 124 Bc + 8A? Be? — 2 A8c3) y$ 
+ (1/9!) (—36A2d— 284 B? + 444?Be — 18A3c?) y® 
+ (1/10!) (1064— 20A Bee + 80A B76 + 20A2 Bc? — 164A2Bed 
+ 64A°B + 76A3c?d — 10A3Be3 — 56.A3c + 2A4c4) y+ ..., (37) 
0,(y) = (1/3!) y+ (1/5!) (2B— Ac) y® — (1/6!) 2Ay® 
+ (1/7!) (836? —3A Be + A*c? + 44d) y’ + (1/8!) (6A2%e — 8AB) y® 
+ (1/9!) (244 Bd — 22A2%cd + 110A? + 4B3 — 64 Bc + 4A? Be? — A3c?) y? 
+ (1/10!) (—80A?d —20A B? + 380A? Be — 12 Ac?) y!? 
+ (1/11!) (160424 5B4— 10A Be + 60A B26 + 10A2B2c?2 + 62428 
— 142A?Bed + Atc* — 52A%c — 5A8Be3 + 6243078) y! 
+ (1/12!) (5804%cd — 72 A? Bc? + 90A2B2c — 8043 — 404 B3 
~ 64A2BS + 20A 4c3) yl? +... (38) 
The solutions must converge for all finite A, B,c and y because of the absence of 
singularities in the differential equation. An inspection of (35), (36), (37) and (38) 
indicates that for the range of the parameters considered, the solutions converge 
rather rapidly. For convenience, a further function ¢,, defined by 
Ps(Y¥) = CP1(Y) + Poly) — OG 3(y), (39) 


is also introduced. It is found convenient in later work to use the set of four 


solutions $5, d3, @, and ¢;. 
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ve The solutions of (18) and (19) are given by linear combinations of the f;(y) 
f and h;(y), these being the values of ¢,;(y) when A and é are assigned the values 
) A,, land A,, —1, respectively. 
Satisfaction of boundary conditions (22) and (26) is insured by writing 
SY) = G2fo(y) +s fo(y) + Oa faly) + 5 f5(y), (40) 
j h(y) = aghg(y) + aghg(y) + ashs(y) +aghq(y). (41) 
> To satisfy the remainder of the boundary conditions, it is necessary that 
») 0 = a,c(A,—A,) +a,ia(4cotd+a?8) —ag, 
| 0 = aef(1)  +a3f3(1)  +agf,(1)+45f5(1), | 
O=ayfi(l)  +agfs(1) +ayfi(1) +a5f5(1), (42) 
0 = a,h,(—1)+a3h,(—1) +a;h;(—1)+agh,4(— 0, 
0 = agh3(— 1) +agh3(—1) +ash;(—1)+agh4( —1)./ 


In order for a non-trivial solution of these algebraic equations to exist, the deter- 
minant of the coefficients must vanish. Therefore we have 





) c(A, — A.) 0 1 ia(4cotA0+a28) —1 
: a, . ' : 
Jx(1) Js) fa) J5(1) 0 
0= — f2(1) fs(1) fa) f5(1) 0 |. (43) 
h(—1) h(-1) 0 hs(—1) h4( —1) 
M-1) Ki-) 0 h(—1) (1) 
‘ This determinant may be expanded into a number of two-by-two determinants 
) by a method such as Laplace’s expansion. This yields the result 


0= A F,+F,A,-H,4%,-hH,+4,F,+ FA; 
+1a(4 cot 6+ a7) (H,F, —F, H,)+¢c(A,—A,) (4, 4;-—H,F;), (44) 


where the following notation has been introduced: 








> fa) fs(1) _| A(—1) 4-1) | 4 
A=|'nq) fay? =| a—1) (-1) 
i) f,(1) H, =| Pl 1) ha(-1) 
fi)’ a | M—-2) —D 
fx) _| 4g(-1)  &(-1) 
fxQ)’ ~* M1): 4=) si 
fx(1) | = h(—1) h3(—1) 
f3(1) e"  t=- 2: 2D 
f5(1) | hg(—1) hs(—1) 
foy | >| a-2 0) 
fs) _ | hg(—1)  hs(—1) 
fs(1) |’ Me =! B(—-1) k(—1) 
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To complete the determination of the eigenvalue c it becomes necessary now 
to approximate the ¢’s, as obviously the infinite series form for the ¢’s must be 
terminated. It was decided to consider a and F to be of the same order, and to 
limit the calculation to fifth powers of « and R (and combinations of them). 
The odd power is introduced by the boundary condition on the normal stresses. 
No special assumption is made concerning c, although it is considered to be of 
the order of one. Computations based on such an approximation prove straight- 
forward but lengthy. To carry out the next-higher approximation (to the sixth 
power) would triple or quadruple the amount of work necessary, which is already 
considerable. 

Treating « and FR to be of the same order gives results which are valid in a 
small circle around the origin in the («, R)-plane. Calculations based on expan- 
sions to the first order in the Reynolds number as made by Yih (1954) would give 
results valid in a narrow strip along the a-axis in the («, 7')-plane. Since the 
results of most stability analyses (including this one, as will be seen) indicate 
the onset of instability with small wave-number, it would seem that an approxi- 
mation valid for a wider range of R# is more desirable. This was the main factor 
in choosing an analysis based on expansion in terms of powers of the co-ordinate y. 

It is desirable to state the range of a and R for which the approximation is 
valid. Certainly if x and R are both much less than one, the results are very good. 
For « and R of the order of one or slightly larger, it is felt that the conclusions 
are still qualitatively true, since the denominators of the terms in the series 
increase in a factorial manner while the numerators increase only gradually. 
However, the prohibitive amount of calculations necessary make any definite 
statements on the range of validity impossible. 

Since the calculations are straightforward, they are not included here. Because 
in several calculations it was necessary to take the difference between two small 
numbers of almost equal magnitude, it was considered useful to leave the 
numbers in fractional form. This also facilitated checking of the results. 

When all calculations are carried out, an algebraic equation in integral powers 
of c is obtained. This is 

0 = —320iz cot 0 +12 (r+ 1) iaR, — 32% cot 6 + $4(r— 1) a?R, cot A 

— 801238 — 712. (r?— 1) PRE + BE (r+ 1) ic3R, + 7142. (r? + 1) ib Ri cot 6 
—5f;ria3 Rj cot 7 + 388(r—1)a4R, cot O— apie 7] 
— S1a5S + Liat R, S(r — 1) 
+ c{ — 1920 —33%aR,(r— 1) — 10240? — *P(r+ 1) a?R, cot d 
+35r0("? + 1) PRY — yGsa2Rir — P(r — 1) ta Ry — 29888 
—35(r? — 1) ia Ri cot 0 — 228 (r+ 1) a4R, cot d— 29 aR, S(r+1)} 
2(256(r + 1) ta R, — 333(r? — 1) o® RE + 192 4ia3 Ry (r + 1) 


+ 8(r2+ oe 2 ia3 Ri r cot O} 


+ cla? Ri[(r? + 1) 182 +297). (46) 








a! 


Vv 
t 


—- &-& ae & ot 


a 
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The root of interest can be computed by assuming ¢ = ¢9+¢,%+¢,R+¢,07+..., 
and substituting this into (46). After this is done, upon separation into real and 
imaginary parts, the result is 
_ #R,(r—1){ 2183 
' 1920 (64,680 
o = %_firt 1) R, 

‘ 1920| 7 


(r+ 1) R, — Sha?S — cot 0 [477 — tt] : (47) 


(12 — 2/3a?) — 16a28(5 — a?) 


af 443,683 — , ane 
(r2+ 1) a? Ri 


~coté 64(5—a8 + ae 
2 ,164,2 


wy 2 POPE . : | 
par — cot? 4[8/3(r + 1) a*Ry]). (48) 


The growth rate o is given by 
= ae,. (49) 
It can be argued (see, for instance, Schubauer & Skramstad 1947) that the 
most likely wave-number to occur is the one with the highest growth rate. This 
value of « can be determined by setting to zero the derivative of o with respect 
to x. The relation obtained by this is 
a (8(r+1)R ; —_ : , 6a! 
ma" | _ ) "(3 — 4a?) — 32Sa?(10 — 3a?) — cot d 64 10 — 4a? + 1 
Ja j \ al, 


443,683 >» 2187 


. er 2PR2 ,; 2 D2 49 32 “ | is 
“mie” (Det rath cotter + TatKy}). (50) 


For the case of vertical inclination, the possible roots of (50) are, besides « = 0, 


a2, = 1[10+ 2 + (100-88 + £2)4], (51) 
where S = =? len 9 (52) 
848 


(In the case of S = 0, the roots reduce to 0, 3.) For small /, (51) may be expanded 
using the binomial theorem to obtain 


a. = 19140-18 + (1— 0-448 — 0-092/2)}, (53) 


) a 


or. for the root of interest, 
ae = P{0-9 + 0-158}. (54) 


5. Discussion of results 
A. The special case of infinite slope 

An inspection of (48) and the corresponding figures 2, 3 and 4 reveals that for 
the case of vertical inclination and zero surface tension, every value of « makes 
c; zero for zero Reynolds number. Hence the «-axis is the neutral stability curve. 
(The analysis does not hold true for large values of «, however, it can be argued 
on physical grounds that in the absence of surface tension there is no restoring 
force when the slope is infinite.) Surface tension has a stabilizing effect (figures 
2, 4) and reduces the range of a for which instability occurs at zero Reynolds 








332 William Paul Graebel 














i ! i 
0 05 1-0 15 2:0 





FicurE 2. Curves of neutral stability for flow between vertical walls with various 
values of the surface tension parameter S. 0 = 0, r= 1-2, 9= 90°. 

















Ficure 3. The rate of amplification of waves of various wave-numbers 
for the vertical case. S = 0, r= 1-2, 9= 90°. 
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Figure 4. The effect of surface tension on growth rate for the vertical case. 
o = 0:25 x 10-3, r= 1-2, 0= 90°. 
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number to the point « = 0, hence it can never completely prevent instability. 
From a physical point of view, the effect of surface tension decreases with 
curvature. From a mathematical point of view, S is always accompanied by a 
curvature term («?), so for very small wave-numbers the stabilizing force would 
be small. 

The shape of the growth rate curves of figure 4 is very similar to the neutral 
stability curve shown by Yih (1954). It might be speculated that the disagree- 
ment between Yih’s and Benjamin’s (1957) results could be due to the difference 
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Ficeure 5. The effect of surface tension on the celerity of the disturbance 
for the vertical case. R = 1-0, r= 1-2, 0= 90°. 


in the approximations which might yield different members of the same family of 
curves. Benjamin’s result would be the true neutral stability curve, while Yih’s 
would be a curve of constant (small) growth rate. 

It is noted from (9) and (47) that the dimensional celerity of the disturbance 
is directly dependent on the difference in densities of the two fluids and the 
velocity of the primary flow. A small Reynolds number implies either large 
viscosity or small difference in densities, hence the speed of propagation of the 
disturbance is small for the approximation considered. The celerity is also small 
for the case of small wave-numbers. The disturbance anticipated by the analysis 
is then in effect almost a standing wave, in the sense that it does not propagate, 
at least not very quickly. 

It is seen from (47) and figure 5 that the surface tension could make the celerity 
negative, hence disturbances would travel uphill. This effect is perhaps un- 
expected, but not unreasonable. Such results have in fact been noted in experi- 
ments conducted by William M. Sangster of the University of Iowa, for layers of 
different thicknesses. 
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If (47) is contrasted with the corresponding results of the stability analyses 
of Yih and Benjamin, it is seen that the celerities represented in their analyses 
have values of the same order as the primary flow for very small Reynolds and 
wave-numbers, while the celerities presented here are many orders smaller than 
the velocity of the primary flow. 

The direct dependence of the disturbance celerity on the difference in densities 
of the two fluids may seem surprising when compared with previous results 
(Lamb 1932, p. 370) for gravity waves. The reason for this difference is that the 
present waves are predominantly influenced by viscosity rather than gravity. 

The fact that the primary flow velocity is zero at the interface may explain 
some aspects of the disturbance celerity. Although a free surface is generally 
regarded as a destabilizing agent, this effect is somewhat diminished in the 
present problem due to the counterflow. If the free surface were replaced by a 
rigid boundary, the problem would be the familiar one of plane Poiseuille flow, 
for which instability is known to occur for Reynolds numbers many orders higher 
than contemplated in the present analysis (see, for example, Lin 1955, p. 28). 

Another interesting feature is that (24) shows that the disturbance amplitude 
is inversely proportional to c. For small c then, even if the y velocity component 
is small, the amplitude of the disturbance may be considerably larger in com- 
parison. 

Naturally, the onset of instability at such low Reynolds numbers does not 
mean the onset of turbulence, but only the onset of waves at the free surface. 
The analysis performed here is, of course, based on the assumption that the two 
fluids do not mix across the interface. 


B. Other slopes 


The results are not as informative for values of # other than 90°. As the inclina- 
tion becomes even slightly less steep, the 4 terms in (48) predominate and exert 
a strong stabilizing influence. For ? more than a degree of two away from ver- 
tical, the results probably are not sufficiently accurate to predict instability, 
and higher-order approximations are necessitated. Neutral stability curves are 
shown in figure 6 for various values of the parameter S for the case @ = 89-5°. 
The curve for S = 0 bends towards the a-axis and predicts instability again 
for all small Reynolds numbers. It is believed, however, that this is not a reliable 
result, and that if higher-order terms were present the curve would be almost 
vertical for small a, and then would go in the positive R, direction, as do the 
curves for non-zero values of S. This same effect is shown in figure 7. The results 
shown in figure 8, with a non-zero value of S, are more reliable in regard to the 
shape of the curve, and show that for inclinations other than the vertical, critical 
Reynolds numbers do exist. It seems reasonable that the critical Reynolds 
number should occur at a = 0. The greater surface tension effect and greater 
dissipation of the disturbance energy at higher wave-numbers would both tend 
to stabilize the flow. Such tendencies were shown for the case of vertical inclina- 
tion. A definite mathematical answer to this matter would require a higher-order 


approximation. 
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Figure 6. Neutral stability curves for the case 0 = 89-5°, for various 
values of the surface tension parameter S. 0 = 0, r= 1-2. 




















FicuRE 7. The rate of amplification of waves of various wave-numbers 
for the case 0 = 89-5°, S=0, r= 1-2. 

















FicurE 8. The rate of amplification of waves of various wave-numbers 
for the case 6 = 89-5°, S = 0-05, r = 1-2. 
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6. Conclusions 

From the results presented, the following conclusions can be drawn: 

(a) For the case of vertical inclination, there are values of « for which the 
flow is unstable at all values of the Reynolds number. The celerity of such 
disturbances is small, but the effect on the interface can be large. 

(6) The stabilizing effect of surface tension is shown, although surface tension 
can never induce complete stability, since it has no effect at zero wave-number. 

(c) For slopes other than the vertical, a critical Reynolds number exists. The 
value of the critical Reynolds number probably is the value of R, which occurs 
at «= 0, although a more thorough investigation is needed in order that a 
definite statement can be made in this regard. 

(d) The presence of the interface is responsible for instability at low Reynolds 
numbers at steep slopes. However, the layer in counterflow does contribute a 
stabilizing influence. 

(e) The critical wave number has been found and is given in equation (54). 
Specific results for the (complex) phase velocity are given in equations (47) 
and (48). 


The author wishes to express his appreciation to Dr C. 8. Yih for suggesting 
the present problem and for his advice during the course of its solution. The 
partial support of this work by the Office of Ordnance Research of the U.S. Army 
is also gratefully acknowledged. 
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Calculations of the motion of non-uniform 
shock waves 


By JAN ROSCISZEWSKI 


Pierce Hall, Harvard University 
(Received 31 March 1959 and in revised form 30 January 1960) 


Problems of the propagation of shock waves and strong detonation waves through 
ducts of variable cross-section and ducts with porous walls, and the interaction 
of a rarefaction wave with a shock and a contact surface in one-dimensional 
unsteady, or two-dimensional steady, flow are solved using a simple rule. Tables 
enclosed in this paper permit efficient calculation of several problems of non- 
uniform shock motion. 


1. Introduction 


The aim of this paper is to give a simple method of analysing various situations 
involving the motion of non-uniform shock waves. The method is a modified 
form of the characteristic rule given by Whitham (1958), or the methods given 
by Chisnell (1957) and the author (1958). The rule developed in this paper enables 
one to treat many cases of the motion of non-uniform shock waves with different 
boundary conditions. Whitham’s characteristic rule is a special case of this 
method when a particular kind of boundary condition obtains. 

The problems treated represent boundary-value problems for a system of 
partial differential equations, in which one boundary condition is given on one 
line, and on another unknown line (the shock locus) a second boundary condition 
is defined by compatibility equations. The rule of this paper enables the behaviour 
of the shock wave to be calculated without solving the system of partial differ- 
ential equations governing the flow. 

Some of the problems discussed below have been solved by the use of methods 
differing from the present one. Some very practicable step-by-step methods for 
difference equations are discussed by Courant & Friedrichs (1948). The problem 
of the propagation of shock waves through a duct of variable cross-section was 
solved by Chisnell (1957), for the case when the gas in front of the shock is at rest. 
Another method of solving this problem was given by the author (1958). In the 
present paper, shock-wave propagation through ducts of variable cross-section, 
when a given non-uniform steady flow exists ahead of the shock wave, is 
discussed. 

The interaction of a simple wave and a shock wave was solved in a closed form 
by Friedrichs (1948) for the case of weak shocks where the entropy changes 
may be neglected. The influence, of a slow perturbation of the otherwise uniform 
piston motion on a shock wave, was discussed by Gundersen (1958). The present 
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method gives good results for arbitrary shock strength and arbitrary piston 
motion. 

A similar procedure applied to the case of plane hypersonic flow around a 
profile gives better results than the shock expansion method.* 


2. General relations for unsteady non-isentropic flow in ducts 


The system of equations governing unsteady non-isentropic flow behind a 
shock wave in a duct can be written in the characteristic form 


Or or (V4 ae a oS aDinA +F(a, V,2,t - 
—=54 ta) = —— + F(a, V, : 2. 
Cat cz 2¢,(k—l1)ea 2 Dt w;t) aad, 
6s = 08 0s a cS aDinA_ . 
— = =+(V-a)— = x3 + Ga, V,2x,t), 22 
ea - Oa ta 
a V a V 
where a ee ee ee 


a denotes the velocity of sound, V the velocity of flow, « and / the curvilinear 
characteristic co-ordinates, S the specific entropy, and A the cross-section area 





— 


x 


FicurE 1. The boundary-value problem for the system of equations (2.1) and (2.% 








bo 
— 


measured along plane, cylindrical or spherical surfaces for the cases of approxi- 
mately cylindrical, wedge-shaped or conical ducts, respectively. 

D @ a 

Dt ct ex 
is the substantial derivative, and F and G are coefficients which express the 
entropy change of an element of the gas and the cross-outflow influence on the 
mass and momentum equations.+ The boundary-value problem for the above 
system of differential equations can be set as follows. In the (x, t)-plane, suppose 
the line LZ, not a characteristic, is given (see figure 1) and that on this line values 


* In an as yet unpublished paper the author has applied similar methods to the cal- 
culation of the shock-wave shape and pressure distribution for hypersonic flow around 
bodies of revolution. 

T See §5. 
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r=r,(t), S=S,(t) and A = A;(x) are known. On the right of a second un- 
known line (the shock wave) the state of the gas is prescribed by V, = V;(2), 
a, = 4,(x) and S, = S,(z). Along this line compatibility conditions must be 
fulfilled. These conditions connect the parameters ahead of the shock wave 
(suffix 1) with parameters on the back of the shock wave (suffix W), and can be 
written in the form 


2 (qm 1 a 
V, =Eil"-y +V,, (2.3) 
v 1/ 
] = 
_ y —_ . 2 2 2 
a a ‘ 
rw =f(W,V,.aq) = ("+ ©) ay (2.5) 
— ([ 2k —. k—1][(k—1) Wi+ 2)"\_ og _ 
Sw = P(W,,S,) = ey In len W I~ E+1}] (k+1)W? IJ S,, (2.6) 
-_ U-Yw as Vy o 
where W, = . Ve = ie etc., (2.7) 


and U is the velocity of propagation of the shock wave. 
We suppose that A = A(x). Then equation (2.1) can be written in the form 


er a es aV olnA : 
on Os a+ ——— = Fla, Va, 6). (2.8) 
Qa 2c,(k—1)ea 2V+a) ca 

Integrating equation (2.8) along «-characteristics (see figure 1) and expressing 
the values on the shock wave by equations (2.5) and (2.6), we obtain* 


rA(W,, Si) a 





Vis — -¢ S 
f( V,,4,)— Tr Js, a(k—I)e,” 
eln Aw ro wy 
dinA = F(a, V,x,t)da, 
JIn 4, 2(J a) J ay, 


where the suffix Z refers to values of the parameters on the L-line. We now apply 
this expression to two neighbouring characteristics 1 and II (see figure 1). On 
subtracting these expressions, we obtain in the limit as characteristic I approaches 
characteristic II,7 


oz i a ( 
i ACses a,)dW,+fr-, (W. py )qV, + fa,( W. Vi.a da, — dr, | 56 ya : 
’ 4? 7h Y 737 a 1 l ; tL IY y ’ 
’ eladinantel wear He, da (9 vS1)-Sz] 
aV aV _if,{ eV , a 
+l rsa. din Ay —- F xv ol, A,+5la| rca) | yg Fy day 


— Fi da, + [dF], (%p—&z), 
where the suffix m refers to the mean value in the interval considered. 
wcor e ° 
* I 3p = = 0, ete., because df = 0 along «-characteristics. 
a 


+ From equation (2.6) we have ¢§ (W,, S,) = 1. 


bo 
bo 
bo 
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For the case when the variations of the flow behind the shock wave are small 
(the linearized case), all values can be written as A = A; +€4,€4=a@,+6,, ete., 
where €, < A,;, €, < a,, etc. We can then omit in the above equation as small 
quantities of the second order all the terms in which the suffix m appears. Then 
we obtain 


fir (W1.V,.a,) dW, +f7-(W1.,4,) dV, + £,(4.4,q) da, 


(dS, + by (W,, S,) dW] + | 


: ; 2(V +a) 


a ‘ 
aS dinA,—F,-da, 
Pre ah |, ‘i iit 


a 
= dy : Fe E= l)c, 


For the general case we may use the formula (2.9) as an approximation. 





-— wan 
S7,,+ A,—F,da,. 2.4 
bs Lt ley cail, dinA,—F,da, (2.9) 


4 


The same procedure can be carried out with the formula analogous to (2.8) 
for the /-characteristics. In this manner we could obtain in addition the value 
of s, = a,/(k—1)—49;, along the L-line. The main difficulty is to find the points 
of intersection of corresponding characteristics with the W-line and the L-line. 
For this purpose an additional assumption is necessary (see § 6 B below). Often, 
instead of r,(t), the value of V;,(t) (for example, the piston velocity) is given along 
the L-curve. In this case, the value of r,(t) depends on s,(t), which can be found 
only if the distribution of r,(t) is known. The situation is now more complicated. 
This kind of boundary condition occurs in the problem of variable piston motion 
behind the shock wave, when the reflected wave reaches the piston. This problem 
was solved by Gundersen (1958), but only in the linearized approximation for 
small variations of the piston velocity when all the flow field can be found in 
closed form.* 

In the special case of boundary conditions such that 7; = const., S; = const., 
A, = Ay = const., we obtain from (2.9) Whitham’s characteristic rule. In the 
case of F = Oand constant state ahead of the shock wave, this rule can be written 
in the form 

fr dU — aw by dU + = | d\n iy = 0. 
+ 4) \qp Ay 
Whitham obtained this rule in a different way: he made use of the equation (2.1), 
in which he replaced derivatives in the characteristic direction by derivatives 
in the shock-wave direction {¢/¢y = ¢/¢t+ U(e/ex)'. He remarks that the slope 
of the «-characteristics is close to that of the shock-wave line (in the linearized 
approximation, they are the same). 

The same result can be obtained (see Chisnell 1957 and Rosciszewski 1958) by 
solving first the linearized case for small shock-wave velocity variations, and 
then, for the case of arbitrarily variable shock-wave velocity, changing the 
parameter in the linearized solution. The quasi-linearized solution obtained in 


* In the general case it is necessary to make use of the expression analogous to (2.9) 
for /-characteristics. Unknown /-characteristics could be replaced by characteristics 
taken from the simple wave solution. In the analogous case discussed in the author’s 
unpublished work dealing with hypersonic flow around bodies of revolution, a similar 
method shows good accuracy. 
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this way is, in fact, the result of applying a series of local linearized solutions to 
the case where finite changes in shock-wave velocity takes place. 

The formula (2.9) will be applied below to many different cases of non-uniform 
shock-wave motion. The results will be compared diagrammatically with finite 
difference calculations obtained from the method of characteristics; these com- 
parisons display the surprisingly good accuracy of the present method. 

A similar boundary-value problem will later be formulated for plane steady 
supersonic flow (ee § 8). 


3. Shock-wave propagation in a duct of non-uniform cross-section 


We shall now apply the rule to the case of shock-wave motion in a tube in 
which the cross-section is initially constant and then varies, and for which there 
is a given constant state ahead of the shock wave (see figure 2). In particular, 
shock-wave propagation in ducts of variable cross-section, when the gas ahead of 
the shock wave is in steady motion or at rest with a variable temperature field, 
can be considered. 

We suppose that A = A(z) and neglect the heat exchange and surface friction 
(F = 0). In the case considered, along the t-axis we have r = const. and A = Ay 
(see figure 2). Taking the f-axis as an L-line, we have dr; = 0,dA, = 0, dS, = 0. 
Now taking a,, V;, S, as variable, and dividing equation (2.9) by dz, we obtain, 
using equation (2.7), 


dU (1- hee (A-0+% a ay ds, 
21 


dx de \1 ~"  N ) dx 2(k—1)c,N dz 
ao roe a ao, (3.1) 
where ul ae ) (3.2) 
R = fir,(W,,V;,, a1) ~ ; et 1+ - , (3.3) 
E = dir(W,,5)) sy i! — : (3.4) 


For the given state ahead of the shock wave, dV,/dx, da,/dx, and dS,/dx are 
known functions of x, and the equation (3.1) is a first-order differential equation 
for U = U(x). In general this equation can be quite easily solved numerically; 
then we obtain the unknown shock wave without having to solve the partial 
differential equations. In the special case of an approximately cylindrical duct 
and small variations of the values of the flow parameters ahead of the shock wave, 
all the coefficients of equation (3.1) can be taken, according to the linearized 
approximation, as constant. Therefore, a closed-form solution can be obtained 
for the shock-wave velocity change and the entire field on the back of the shock 
wave. The discussion of this solution was given in the author’s previous work 
(Rogeiszewski 1958). 
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In the more general case when the time-dependent field ahead of the shock 
wave is given [V, = Vi(z,t), a, = a,(x,t), S,; = S,(x,t)] and, as before, dr, = 0, 
dA, = 0, dS; = 0, in place of d/dx on the right-hand side of equation (3.1) we 


must put ¢/cx+(1/U)(c/ct). We then obtain the equation for the shock-wave 
velocity in the form 


dU * 
dx 


F(U, x,t 


~— 


~ F(U,2, [’ = 
~ 0 





This equation can also easily be solved numerically. 
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FIGuRE 2. Propagation.of a shock wave through a duct of variable cross-section, (a) sub- 
sonic flow behind the shock wave, (b) supersonic flow behind the shock wave. 


The solution of the problem when time-dependent boundary conditions on 
curve L are given [r,; = r,(t), A; = A,(t)] is also possible, and it is not very 
complicated when the curve L lies outside the domain of dependence bounded 


by the /#-characteristic and the shock wave, i.e. when the condition on the 
L-curve does not depend on the shock-wave motion. 
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In the special case when a, = a) = const., V, = 0, S; = const., using equations 
(2.3)-(2.6) we obtain from equation (3.1), 


A rU | 2(k+1) k+1 | 

In-=-=- = = += - 

Ay Jo, \V{2(k—1) (kU? —(1/U2)] + 6k -k?-}} U -(1/U)) 
x/p_-2w* lao, (3.5) 


\(" 2(k-1)J 


where R and E£ are functions of U given by equations (3.3) and (3.4). 


rU rU U 
0 -| {3dU0 U ~f £348 0 -| {}d0 
2 J2 2 
1-2 —4:514 3-5 3-186 5:8 5-817 
1:3 — 3-563 3-6 3-338 5-9 5-854 
1-4 — 2-612 3-7 3-481 6-0 5-991 
1-5 — 1-999 3-8 3-623 6-1 6-076 
1-6 — 1-490 3°9 3-757 6-2 6-160 
1:7 — 1-053 4-0 3-891 6-3 6-242 
1-8 — 0-668 4-1 4-017 6-4 6-323 
1-9 — 0-322 4:2 4-144 6-5 6-402 
2-0 0-000 4:3 4-265 6-6 6-482 
2-1 0-294 4-4 4-385 6-7 6-559 
2-2 0-587 4:5 4-501 6:8 6-635 
2°3 0-845 4-6 4-616 6-9 6-710 
2-4 1-102 4:7 4:727 7-0 6-785 
2°5 1-330 4:8 4:838 7-1 6-858 
2-6 1-559 4-9 4-944 7-2 6-930 
2-7 1-767 5-0 5-050 7-3 7-001 
2-8 1-974 5-1 5-152 7-4 7-071 
2-9 2-160 5:2 5-254 75 7-140 
3-0 2-356 5-3 5-352 7-6 7-208 
3-1 2-537 5-4 5-449 77 7-275 
3°2 2-707 5°5 5-543 7:8 7-341 
3-3 2-870 5-6 5-636 7-9 7-407 
3-4 3-033 5-7 5-727 8-0 7-472 


TABLE l| 


In the limiting case when U -> 0 we obtain the asymptotic formula 


2\, U 
= (faye + 5) ng. (3.6) 


In fact, equation (3.5) is the Chisnell answer presented in another form. Some 
numerical differences between Chisnell’s tables and present formula are, however, 
found (see the converging cylindrical wave calculation). Values of the integral 
in (3.5) are given in table 1 for k = 1-4. 

Equation (3.4) enables us to examine the influence of entropy changes on 
shock-wave motion. The coefficient depends on 


1 d(Sy—S;) 
c, dW, 


p 


E 
= 45 
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where F given by equation (3.4) expresses the entropy: change influence. When 
W, = 1 (very weak shock wave) and W, = 00 (very strong shock wave) E vanishes, 
This means that entropy changes for hypersonic shock motion can be neglected. 
Table 2 gives values of £ as a function of W, for k = 1-4. 

The results of calculations by the use of Table 1 are compared with step-by- 
step calculations* applied to the basic differential equations for various duct 
shapes and various initial shock strengths (figures 3-10). In special cases these 


W 0 1-5 2 3 4 5 6 8 10 15 20 00 
E 0O 0-097 0-185 0-246 0-241 0-221 0-198 0-161 0-123 0-092 0-0714 0 


TABLE 2 
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Ficure 3. Comparison of the present method with step-by-step calculations by the 
method of characteristics for the case of shock-wave propagation through a duct of vary- 
ing cross-section area given by A/A, = e%* 0) (divergent duct) and initial shock strengths 
(a) Uo/a) = 5 and (6) U,/a) = 3. The shock is assumed plane. —-—-, Step-by-step calcula- 
tion; ———, present method. 


* This calculation was based on the method of characteristics. The accuracy was checked 
by taking some mean characteristics which give the errors lying within the limits of 
accuracy of the wave diagram calculations. 
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results are compared (figure 9) with Payne’s (1957) calculations for a cylindrical 
| convergent wave. In this case, because of other boundary conditions (in Payne’s 
work V(1,t) + const.) Payne’s results differ from those given by the present 
method and from those given by finite difference calculations by the method of 
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Figure 4. Comparison of the present method with calculations by the method of charac- 
teristics for the case of shock-wave propagation through a duct varying cross-section area 
given by A/A, = e?(?/*) (divergent duct) and initial shock strength U,/a) = 7. The com- 
parison is made for both a plane and spherical shock wave, originating at «= 0, for 
» present method. 





x/x) > 0-1. —-—-, Step-by-step calculations; 
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Figure 5. Comparison of the present method with calculations by the method of charac- 
teristics for the case of shock-wave propagation through a duct with A/A, = e) (diver- 
gent duct) and initial shock strength U,/ay= 1-5. —-——-, Step-by-step calculations; 
, present theory. 
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characteristics. The results are not strictly insensitive to the initial conditions, 
see Whitham (1958). Chisnell’s tables show some disagreement with the present 
tables, and they give results very close to Payne’s calculations for different 
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FicurE 6. Comparison of the present method with calculations by the method of charac- 
. , — S. ee 
teristics for the case of shock wave propagation through a duct with =exp| —{—+=,4 
Ay Ly 229 
(convergent duct) and initial shock strength U,)/a) = 2-5. —-——, Step-by-step calculations; 
—, present method. 
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Ficure 7. Comparison of the present method with calculations by the method of charac- 
fs : ow al a xy? 
teristics for the case of shock-wave propagation through a duct with — = exp] — +- ) 
Ay Xo 2x9 
(convergent duct) and initial shock strength U,/a,= 1-3. O, Step-by-step calculations; 
x, present theory. 


boundary conditions. For other cases calculated in this paper, Chisnell’s tables 
also show some deviations. Results of the comparison of the present method with 
experimental results and calculations by Hertzberg & Kantrowitz (1950) are 
given in figure 10. 
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4, Propagation of strong detonation wave through ducts of non- 
uniform cross-section 


Here we shall be concerned with the propagation of a strong detonation wave, 
i.e. the detonation wave described by the points on the Hugoniot curve lying 
above the Chapman—Jouguet point (see Courant & Friedrichs 1948). Such a 
detonation wave can be obtained by additional compression behind the detona- 
tion wave, for example, by piston motion or by the propagation of the Chapman-— 
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FicureE 8. Comparison of the present method with calculations by the method of charac- 

teristics for the case of a very strong shock wave (U,/a,) = 20) propagating through a duct 

of variable cross-section A/A, = e**/*), The gas is treated as ideal with k= 1-4. @, Step- 

by-step calculations; O, present method. 


Jouguet detonation wave through a duct of suddenly convergent cross-section. 
We assume the strong detonation wave to be propagated through a duct of 
initially constant cross-section, which at some point changes into a duct with 
variable cross-section. The gas ahead of the detonation wave is assumed to be 
at rest with a constant velocity of sound. 

In a divergent duct the velocity of a detonation wave can be decreased only 
to that corresponding to the Chapman—Jouguet wave. The only difference 
between the solution of the previous case and the present one is in the different 
compatibility equations for shock and detonation waves. 

We assume different values of k = c,/c, in front (suffix 1) and at the back 
(suffix W) of the detonation wave. 
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FIGURE 9. Comparison of the present method with calculations by the method of charac- 
teristics (assuming V(1, ¢) = const.) and Payne’s calculations with different boundary 





conditions (V(1, ¢) + const.). —--—-—, Step-by-step calculations with boundary condition 
V(1, t) =const.; ——-—-, finite difference calculations by Payne V(1, t) + const.; - 
present method ; , Chisnell (published in Payne’s work). 
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Figure 10. Comparison of the present method with Hertzberg & Kantrowitz’s experi- 
mental curve and their calculations by the method of characteristics. O, Experimental 
points (Hertzberg—Kantrowitz); ———-, finite difference calculations (Hertzberg—Kantro- 
witz); ———, present method. 
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Sed 


The equations of conservation of mass, momentum and energy (see Courant 


& Friedrichs 1948) giv 
- . 1 - ] l 2(k,—-1)Q 2(k, — k,) 
if = mit 2 : 1 a. 
~ hy wai *y)+ \(2 k ine +1) * (ky +1) a ky(ky — 1) (ke + 1)J 
(4.1) 
where U = U/a,and V = Vjz-/a); Q is the heat released by the combustion process 
Also, 
ee ] 
9 
a= |\eo- V) (P+i a) (4.2) 
ae —Ina— i: In a =, (4.3) 


where @ = ayy/dy 
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Present method 


Comparison of the present method with calculations by the method of charac 
teristics for the case of propagation of a strong detonation wave in a duct of cross-section 


FIGURE 11. 
given by A/A, = e#?'20), with k, = 1-4, k, = 1:3 and Q/ag = 28-3 


step by-step calculations 
From equation (2.9) (dV, = da, = dS, = 0) we obtain 
% A e "3 V +a) (1 qV l ll a 
Ag J Us Va \2dl ky—1 ol cVadl 
siege sla) | at oe |lav. (4.4) 
cU e@VdU ca\eU @eVadU!)) 


where the derivatives in equation (4.4) can be calculated by the use of equations 


2) and (4.3). 


(4.1), (4.2 
Comparison of the calculations by the use of the present method and step- 
by-step calculations are given in figure 11 for k, = 1-4, k, = 1:3 and Q/a3 = 28-9 
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5. Shock-wave propagation through a duct with porous walls 

In this section we consider the shock-wave attenuation due to the cross- 
outflow through porous walls of the duct. We take a cylindrical duct and we 
assume the cross mass flow 7d,/pv to be small in comparison with pV A, where v 
is the cross-flow velocity, 7 is the ratio of the area of the pores to the area of 
the wall, and ¢ = vy. 


(a) ie ae 


S= const. Soi 2: A 
eo, ~~ Shock wave _ 
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FicureE 12. Flow induced by shock-wave propagation through a duct with 
A I 
porous walls (7<1). (a) The case V°< a; (b) the case V°>a°. 





We apply the basic mechanical equations to the ideal flow through the duct 
with the porous wall (figure 12). We assume all parameters to be constant in the 
given cross-section. The continuity equation is 


. Parr o ft 
Pi Vi Ao-pPr Vit Ao —7d, pedx = ay pA,dz, 
7 Ty] c o ry 
or for x; > 2, Dp ay app 
+ ( re —_ e 5.1 
Dt * ex d, mm 


For the momentum equation, we assume that the cross-stream has no momen- 
tum in the direction of the axis of the tube. Then we have 


ar 
C 


(Pr- Py) Ag+ py ViAo—Pr Vir Ao = a 


pV A,dz, 


t Ja, 
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or, using equation (5.1) for x, > x,;, we obtain 


lM sa (5.2) 
Dt pcx dy 


The energy equation is 

_ (Vik p V2 k p 

| I ee 1 Pu _ 
Py ui 2 tea 0 —pu¥al-5 + ET pn) Ae 


é 7 si : oo, 
5 aa €. teh 


or, using equations (5.1) and (5.2) for x; > a;, we get 








1 DS _ 20 
c, Dt ad 


|=), M2 “Sel. (5.3) 
Pp 

We shall calculate the velocity of the cross-flow by assuming quasi-steady flow. 
This is based on the fact that due to a sudden change in the cross-section the local 
derivatives (in the equation of motion), are small in comparison with the con- 
vective terms. Then we apply the Bernoulli equation and assume the cross-flow 


to be isentropic 


_¢ || 2k p Pe) - ~ | [| 2a? Pa . 4 5 
. =o Hmibl “—) é/ea|!- (7) } —s 


where p, = const. is the outside pressure, p = p(x,t) is the pressure in the duct, 
and € is the loss coefficient. 
The pressure in this formula can be calculated from 


Pp a\?kk-)  — (S3-S8 ici 
a= I> expt ——_— Ff. (5.5) 
Pe as Cy aaa C, 


Taking into account the three equations of motion (5.1), (5. 
tion (5.4) we find the coefficients in equations (2.1) and (2.2) 
entropy change and the cross-mass flow influence 


2), (5.3) and equa- 
which express the 





» - p Pt  — p a 
a — |—? 1—-M+ M?+ 7) 1 , 
. © Neal} =I Hl 2 ai o) ' 
(5.6 
2a2nt /{ 2 p el a ae Dp, \* ") 
G=-—- : 1—(-" 1+M+ M?2+¢2}1-(44 , 
dy f rer “1 | . 2 | (") 


where M = V/a. 
We obtain from (2.9), where we substitute dx, = 0 (the t-axis is taken as the 


L-curve) and 
dx 


Vp + Oy 


[ at = [nae, (5.8) 


day = day = 
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ry 





where BU) = 7 2 aE ] 
( w+ ay) | 2(k —1) 


E and R are given by (3.3) and (3.4), respectively. The integral of the left-hand 
side of equation (5.8) is given in table 3 for k = 1.4 and € = 1. 

In figure 13 a comparison of the present method with the numerical method 
of characteristics is given for 7 = 0-1 and U, = 4. In this case the numerical 
method was not very accurate. 


oa po ag a rU dU asf U dU 
l l = 

Jo, BU) Jom La ) 
1-2 — 2-775 3:6 0-755 6-0 1-295 
1-4 — 0-900 3:8 0-814 6-2 1-328 
1-6 — 0-425 4-0 0-870 6-4 1-360 
1-8 — 0-172 4:2 0-922 6-6 1-391 
2-0 0-000 4:4 0-972 6:8 1:42] 
22 0-144 4-6 1-019 7-0 1-450 
2:4 0-264 4:8 1-064 7-2 1-478 
2-6 0-367 5-0 1-107 7-4 1-506 
2-8 0-460 5-2 1-148 7-6 1-532 
3-0 0-544 5-4 1-187 78 1-558 
3:2 0-621 5:6 1-224 8-0 1-584 
3:4 0-691 5:8 1-260 


TABLE 3 
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FicurE 13. Comparison of the present method with the calculation by the method of 
characteristics for the case of shock-wave propagation through a duct with porous walls. 
The ratio of pore area to wall area 7 = 0-1. Initial shock strength U)/a) = 4. ©, Step-by- 
step calculations; O, present method. 
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6. Interaction of a rarefaction wave and a shock wave in one-dimen- 
sional unsteady flow 
A. Case of head-on collision 

Until now, these problems for arbitrary waves were solved by the use of finite 
difference methods. A discussion of these methods is given by Geiringer (1948).* 
We discuss the case of interaction of a simple wave and a shock wave of arbitrary 
strength (figure 14) in a duct of constant cross-section.t The flow behind the 
shock wave is governed by the unsteady non-isentropic flow equations. 


‘| o 






a -characteristic 





— 
a x 





FicurE 14. Head-on collision of a shock wave and simple rarefaction wave. 


Introducing a simple wave relation ahead of the shock wave 
da, = —3(k—1)d¥,, (6.1) 
and taking into account the fact that «-characteristics cross the uniform flow 


region behind the shock wave (see figure 13) (dr; = 0), we get from equations 
(2.9) the first-order linear ordinary differential equation 


. Ir ip 1 = "—a Tr ~~ 
! _ x) + Ah 7] ae = ¥(M), (6.2) 
1+4(k-1)V, JdwW, 
where from equations (2.4)—(2.7) 
— ae = 
x(W)) eT k+l! (W,- 7) tan (6.3) 
4(47, +=) 

u(W,) = —"_ EB (1 + ae ee (6.4) 
7 ( 1 Be l 4 k+ l Ww (k+ 1}? a, . 


where a,, and £ are given by (2.4) and (3.4) as a function of W,. Table 4 contains 
the values of y(W,) and y(W,) for k = 1-4. 


* The numerical example in this paper is incorrect. The principal equations of motion 
are not satisfied for the shock wave before interaction. 
+ This case is a special case of §3 when A = const. and V, = V,(z, t), ete., but we solve 
this here in a simpler manner. 
23 Fluid Mech. 8 
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The results for the special case of V, = 0 and W¢ = 4 are compared in figure 15 
with a step-by-step method of calculation (method of characteristics) applied to 
the basic non-linear partial differential equations. In this case the accuracy of 
the present ‘method is excellent. The numerical calculations are very difficult 
because of the slow convergence of the finite difference method. Five iterations 











Wy x(W,) y(W,) Wy x(W,) v(W,) 
1-0 1-00 3°32 4:6 2-92 1-58 
1-2 1-11 2-82 4-8 3-10 1-58 
1-4 1-21 2-50 0 3°22 1-57 
1-6 1-32 2-28 5:2 3°33 1-56 
1:8 1-43 2°12 5:4 3°44 1-55 
2-0 1-54 2-00 5:6 3°56 1-54 
2-2 1-65 1-90 5:8 3°68 1-54 
2:4 1-76 1-84 6-0 3°80 1-53 
2-6 1-86 1:77 6-2 3°92 1-53 
2-8 1-98 1-72 6-4 4:03 1-52 
3-0 2-08 1-69 6-6 4:15 1-52 
3-2 2-19 1-66 6:8 4:27 1-51 
3°4 2-30 1-64 7-0 4:38 1-51 
3°6 2°41 1-63 7-2 4-51 1-50 
3:8 2-52 1-62 7-4 4:63 1-50 
4:0 2-64 1-61 7-6 4-74 1-49 
4-2 2-76 1-60 7:8 4°86 1-49 
4-4 2-88 1-59 8-0 4-98 1-48 
TABLE 4 
4-8 
e 
1-6 
ra | 
1-0 l L | ] 
0 0-2 0-4 0-6 0-8 


V,/a, 
Figure 15. Comparison of the computation by use of the present method with the com- 
putation by the method of characteristics for the case of head-on collision of a shock wave 
of initial strength W$/a$ = 4 anda simple rarefaction wave. @, Step-by-step calculations ; 
O, present method. 
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were made in each step for obtaining points of intersection of the characteristics 
a and the shock wave in the wave diagram; calculations in this case were then 
accurate. 
B. Case of the merging of a rarefaction and shock wave 

We consider the interaction of a rarefaction wave and a shock wave propagating 
in a duct of constant cross-section with a constant state ahead of the shock. 
Behind the shock wave is a piston which initially has a constant velocity (see 
figure 16). |: 
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\ /S# const. 


Piston 
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Simple wave 











= _— 
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FIGURE 16. Merging of a shock wave and a rarefaction wave 
caused by non-uniform piston motion. 


In the present case we have variable boundary conditions given on the 
arbitrary L-curve lying near the piston path. We assume that no reflexion wave 
going back from the shock wave reaches this curve. 

Along the given curve we have the simple wave relation 


a V, 
a ‘i - t = const.; 
therefore dr, = dV, = dV, (6.5) 


where V,, is the variable velocity in the simple wave prescribed by the piston 
motion. 
Using equation (2.9) we obtain 


a0|R- a | 


a(k—1) —4y(U)dU = dv, (6.6) 


Ill 


where R and E are the functions of U given by equations (3.3) and (3.4). 
After integrating we obtain the following equation 


Table 5 contains values of the integral in equation (6.7) for k = 1-4. 


23-2 
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In this case the entropy change plays an essential role; and the shock expansion 
method does not give a good approximation. Equation (6.7) only enables the 
calculation of a change in the shock-wave velocity caused by the piston-velocity 
variation. To obtain the history of the shock-wave motion it is necessary to find 
the shape of the characteristics «. 


1f¥ zs ‘ . 1(U ae os 1U ee 
TC 5] ywo)a O Al w(T) dl C = y(T) a0 
«JU, “JU, “JU, 

1-2 — 0-949 3-2 1-190 5:2 2-795 
1-3 — 0-806 3° 1-274 5:3 2°872 
1-4 — 0-663 3:4 1-358 5:4 2-950 
1-5 — 0-543 3°5 1-44] 5-5 3:027 
1-6 — 0-424 3-6 1-523 5-6 3°104 
1-7 — 0-314 3°7 1-605 5:7 3°181 
1-8 — 0-205 3°8 1-686 5:8 3°258 
1-9 — 0-102 3-9 1-767 5-9 3°335 
2-0 0-000 4-0 1-848 6:0 3:412 
2-1 0-098 4-1 1-928 6-1 3°488 
2-2 0-196 4-2 2-008 6-2 3°565 
2-3 0-290 4:3 2-088 6:3 3°641 
2-4 0-385 4-4 2-167 6-4 3°718 
2-5 0:477 4:5 2-246 6°5 3°794 
2-6 0-568 4:6 2°325 6-6 3°870 
2:7 0-707 4:7 2-404 6-7 3:947 
2°8 0-847 4:8 2-483 6:8 4-023 
2-9 0-923 4-9 2-561 6-9 4-069 
3-0 1-020 5-0 2-639 7-0 4-115 
31 1-105 51 2°717 


TABLE 5 


On the basis of the step-by-step calculations it is observed that the «-charac- 
teristics are not very different from straight-lines. 

Using these results we can calculate the change of the shock strength along 
the tube axis. Comparisons of the calculations by the present method and the 
finite difference method are given in figures 17-19 for various initial shock 
strengths and various modes of piston motion. 


7. Interaction of a rarefaction wave and a contact surface in one- 
dimensional unsteady flow 

The present method enables us to obtain in closed form the solution for the 
interaction of a rarefaction wave and a contact surface in an unsteady one- 


dimensional flow. In a similar way, the analogous problem in plane steady flow 
can be solved. Until now, the closed form solution only enabled the calculation 
of the finite value of the contact surface velocity (after interaction) to be done 
(see Billington 1956). The full history of the contact surface motion was obtained 
by step-by-step calculations (see Courant & Friedrichs 1948). 
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We assume the flow to be isentropic on the two sides of the contact surface 
(figure 20). The compatibility equations on the contact surface are 
Vie = Voc = Ves) ~ 
f (7.1) 
Pic = Poe = Pe) 


(no velocity and pressure jump take place on the contact surface). 
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(b) 
FicurE 17. Comparison of the present method and calculations by the method of charac- 
teristics for the case of a shock wave of initial strength U,/a, = 4 merging with a rare- 
faction wave caused by a suddenly stopped piston. (a) Shock-wave velocity versus 
piston velocity. (6b) Change of the shock-wave propagation velocity, straight line e- 
characteristics being assumed. —-——, Step-by-step method ; ——, present method. 
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From the second relation (7.1) and the equation of isentropic flow, we get 


2k, (k,—1) 2k,/(k,—1) 
(=x) ai =. (>) = (7 2) 
a} ay, ‘ : 


where a!, af are the initial values of the velocity of sound on the right and left 
side of the contact surface. 
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Figure 18. Comparison of the present method and calculations by the method of charac- 
teristics for the case of a shock wave of initial strength U,/a) = 3 merging with a rare- 
faction wave caused by suddenly stopped piston. (a) Shock velocity versus piston velocity. 
(b) Change of the shock-wave propagation velocity, straight line a-characteristics being 
assumed. ———-, Step-by-step method; ——, present method. 
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We have ds, = 0 for the transmitted wave (simple wave) and therefore along 


the contact surface 
a, _ YX, 


ae Ve 


—08 bam 


const. 
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FiGuRE 19. Comparison of the present method and calculations by the method of charac- 
teristics for the case of a shock wave of initial strength U,/a, = 4 merging with a rare- 


faction wave caused by the retarded piston motion dV,/dt =W, = —30(a /ty). (a) Shock 
velocity versus piston velocity. (b) Change of the shock-wave propagation velocity, 
straight a-characteristics being assumed. -~—-—-, Step-by-step method; ———, present 


method. 
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For this case, in an analogous way to that leading to equation (2.9), we can write 

the following equation 

dV, dag, 
ae ee 


dr, = = dr, = dV, (7.4) 


Using equations (7.2) and (7.3) and integrating, we obtain the following 


expression ,, ,, - : ve : i 
I V.-ve ao [V,-V? (K_—1)(ky—1) (Key /k) _™ = 
— (ky = V,— V5. (7.5) 


2 kg—1| 2a 
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FicurE 20. Time history of the interaction of a rarefaction wave with an air-helium 

contact surface, computed by the present method with the assumption of straight 

a-characteristics and by the method of characteristics. , Step-by-step calculations ; 
, present method. 


This equation enables one to calculate the velocity of the contact surface for 
a given change in the simple wave velocity V,, = V,(a/t). Assuming, as before in 
the case of the interaction of a rarefaction wave and a shock wave (see § 6B), 
that the «-characteristics are straight lines, we get the history of the contact 
surface motion. The results of the calculation are compared with the results of 
the method of characteristics and presented in figure 20 for an air-helium contact 
surface. 


8. Wave interactions in plane steady supersonic flow. The basic 
expression for plane flow 


In a way similar to that used for unsteady one-dimensional flow, we introduce 
the analogous expression for plane flow. We make use of the equation in the 
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characteristic direction by assuming that the curvature of the characteristics 
is small (see Howarth 1953, p. 75): 


a 


f= 0. (8.1) 


OV, U,OUg  v 
C lo 


. C 
) ie Oo ees 
H Uzg0H purge 
Here, « and f are the characteristic and the orthogonal curvilinear co-ordinates, 
respectively, p is the pressure, and v,, v, are the velocity components in the 
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Figure 21. The boundary-value problem in plane steady supersonic flow. 


and # directions, respectively. We take the ~-co-ordinate along the C., charac- 
teristic. The compatibility. equations for an oblique shock wave can be written 
in the form (see Stanjukovitch 1955) 
, , cos D4 
eat . 8.2 
WW loos (y—0) ( ) 


where y, @ are defined in figure 22, 


Dy 2k eee k-1] "er 
2 ies k—1 tan (y—9) 
= M¢sin?} : me 8.4 
oe 1 51n iE tani y ( ) 


From these equations we obtain 
‘ , COSY COS fly : . 
Uw = Vy C08 se = | . -=f(y,V,), 8.5 
alt WW 2 1 cos (y —9) fQ 1) ( ) 
Uaw = Vy Sin lg = P(y,V), (8.6) 
Pw = (7%), (m7) 


where the Mach angle “, = sin~!1/M,. We assume that we can express all the 
parameters of the flow in terms of the velocity V, in front of the shock wave. The 
boundary-value problem can then be presented in a similar way as in §2. We 
assume that we are given an arbitrary curve L which is not a characteristic (it 
can be the surface of a body, see figure 21), and that on this curve the boundary 
conditions V = V,(x), p = p,(2) are given. 
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As before in §2, we integrate the expression (8.1) along a C., characteristic 
between the L-line and the shock wave. Expressing the values on the shock wave 
by equations (8.5), (8.6) and (8.7), we obtain 


PAY, Vi) a ry(y, Vi) v, 
P ~ 7 . 2 ‘ — 
tiv A De Me = dug + dp = 0. 
“UBL Up ‘Pr PY 


Applying the above equation to two neighbouring characteristics C _,; and C1, 
(see figure 21), and subtracting one of these relations from the other, we obtain 
in the limit as C,; > C.1 
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FIGURE 22. Head-on collision of an oblique shock wave and 
a Prandtl—Meyer expansion in steady plane flow. 


Here values with the suffix m are the mean values in the interval of integration. 
In the linearized case, all the parameters can be written as V = V;, +e,, etc., 
where €, < };, etc. We see that the terms with mean values are small quantities 
of the second order. Neglecting these, we obtain 
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We shall henceforth apply equation (8.9) as an approximation for cases in which 
linearization is not valid. This equation is analogous to equation (2.9) for 
unsteady flow. 

We shall now make use of this expression to solve the problems of the inter- 
action of a Prandtl—-Meyer expansion and an oblique shock wave, and plane 
hypersonic flow around a body. 


9. Head-on collision of an oblique shock wave and a Prandtl-Meyer 
expansion 

In a similar way to that for the unsteady wave interaction discussed in § 6, we 
now apply the relation (8.9) to the head-on collision of a simple rarefaction wave 
(Prandtl—-Meyer flow) and an oblique shock wave (figure 22). 

In front of the shock wave, the following equations are satisfied by the Prandtl-— 








Meyer flow ", ae z az 9 1) 
2 k-1 k-V’ i“ 
Pr = ae (9.2) 
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With these relatiors, we obtain from equation (8.9) the following expressions: 
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Using the compatibility equations for the shock wave (8.2), (8.3) and (8.4), 
and the equations (9.1) and (9.2), we can express all coefficients of equation (9.4) 
as functions of y and V, alone. We then obtain a first-order ordinary differential 
equation for y as a function of V,. Numerical calculations on the basis of equation 
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FiGuRE 23. Comparison of the present method with the method of characteristics for the 
case of head-on collision of an oblique shock and a Prandtl—Meyer expansion. 


(8.4) enable the computation of the shock-wave path without solving the partial 
differential equations governing the flow on the back of the shock wave. Equation 
(9.4) could be used to construct suitable tables. 

A comparison of the results obtained with those of finite difference calculations 


« 


is presented in figure 23 for M? = 3, y® = 52°. 


10. The shock wave in plane flow around a body 

We consider the problem of the calculation of the shock-wave shape in plane 
supersonic flow around a body (see figure 24). We assume the shock wave to be 
attached (i.e. the nose angle of the body is smaller than the critical one). This 
problem corresponds to the merging shock and rarefaction wave in one-dimen- 
sional unsteady flow considered in § 6B. 

Taking the body surface as a curve L and assuming that there is no reflected 
wave from the shock wave (a small reflexion is observed according to Egger’s 
et al. (1955) paper), we can calculate all the parameters on the body surface from 
the Prandtl—Meyer flow relations. From (8.9) we obtain 
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where B and P, are given by equations (9.5). 
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Expressing all the parameters in the denominator of equation (10.1) in terms 
of y and in the numerator in terms of 6, we get the following equation 


Y poz 


Dy) dy = 


0 
“7 J %, 


(0) dO, (10.2) 
where ®(y), U(@,) are functions given by the denominator and numerator in 
equation (10.1), respectively. 

The equation gives the variation of the shock-wave angle y as a function 
the body shape given by 0,, taken on the same characteristic. 
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FiGURE 24. Supersonic plane flow around a body (the calculations by the method of 
characteristics and the shock expansion method are taken from Eggers et al. 1955). 
—, Present method ; ———— , method of characteristics ; —-—--—, shock-expansion method. 


This method differs from the shock-expansion method, where all the para- 
meters are assumed to be constant along the C, characteristics (i.e. a simple 
wave). In the present solution @ varies along the characteristics. These solutions 
are the same only in the linearized case. However, the assumption that the C., 
characteristics are straight lines was made in the calculation of the shock-wave 
shape. 

Figure 24 presents a comparison of the present theory, applied to a 10% 
thick biconvex airfoil at 1, = «©, with the results published in Egger’s et al. paper 
obtained by the method of characteristics. (This comparison is not very accurate 
because the figure drawn in their paper is not precise.) 

An interesting problem is the influence of the entropy change on the flow 
behind the shock wave. We obtain for oblique shock waves (in a similar way as 
in §3) the following formula 


1 dS ( M, sin y l | ’ 
-= 4 ‘ - M,cosy. (10.3 
c, dy \2kM?sin?'y—k+1 Msiny[(k—1) M?sin?y + 2] eer 
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For M,= a 
tds 2 ai 10.4 
-_= ot y. : 
c, dy k Y ( ) 


We then get for a given curvature of the shock wave a weak entropy gradient 
behind the shock wave, and hence low vorticity for y close to 47 (a detached 
shock wave). 

For the given profile, it follows from the above method of shock-shape cal- 
culation that the higher the Mach number ,, the smaller is the curvature of 
the shock wave. It therefore follows that the entropy change at higher Mach 
numbers is smaller because of the smaller curvature of the shock wave. 


11. Discussion 

The results obtained in the present paper seem to be very useful for many cases 
of non-uniform motion of shocks, detonation waves, and contact surfaces, but 
the present method enables us to find only the position of the unknown dis- 
continuity line. In the case of small disturbances, when the linearized approxi- 
mation can be applied, it is possible to find the complete flow field. In the general 
case, it is necessary for the determination of the flow field to solve the modified 
Cauchy problem for the basic partial differential equations with the boundary 
conditions given on the known discontinuity line. 

The results obtained in this paper are very close to those obtained by step-by- 
step calculations based on the method of characteristic, because the neglected 
expressions containing mean values in equations (2.9) and (8.9) have small 
contributions. This is because the curvature of the «-characteristics is small 
in all the problems considered. The assumption of straight characteristics, which 
was made in some of the cases, gives results which differ little for both methods. 
This difference depends on the entropy change behind the shock wave. 

In this paper the results obtained were compared with step-by-step calcula- 
tions mainly for shock strengths U = 3, 4, 5, where the influence of entropy 
changes is greatest. For weak and strong shocks the results are therefore better 
than for shocks of medium strength. 


The author wishes to express his gratitude to Mr P. Kijkowski and Mr §. 
Pietrzyk for their help in carrying out the computations necessary for the ex- 
amination of the present method; and also to the Fluid Mechanics Department of 
the Institute of Basic Problems of the Polish Academy of Sciences for financing 
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Velocity and enthalpy distributions in the compressible 
turbulent boundary layer on a flat plate 


By D. A. SPENCE 


Royal Aircraft Establishment, Farnborought 


(Received 1 November 1959) 


The object of this paper is to present a unified account of the distributions of 


velocity, shear stress and enthalpy in the compressible turbulent boundary layer 


on a flat plate. As a start, the set of velocity profiles measured over a range of 


heat-transfer conditions at Mach numbers between 5 and 8 by Lobb, Winkler & 
Persh (1955) is examined. It is found that by plotting in terms of the Howarth 


. 


y 

variable 7 = | (p/p,)dy, the outer parts of the profiles for different Mach 
7 0 

numbers are brought together on a single curve of the approximate form 


u/u,, = (n/A)'", A being the transformed boundary-layer thickness. By evalu- 
ating the reference density p, and kinematic viscosity v, at the so-called ‘inter- 
mediate enthalpy (Eckert 1955) the inner parts of the profiles can also be collapsed, 
although less completely, to fit a ‘law of the wall’ w/w, = A log (yu,/vy)—c) +B. 
Here w, = (7,,/))?,and A, B and care the same constants asin incompressible flow. 

These properties provide a physical starting point from which the remaining 
features of the mean flow can be calculated. By substitution of appropriate 
stream functions in the equation of motion the distribution of shear stress 7 in 
inner and outer regions is found; this approximates to the form 

T/T, = 1—(u/u,.)™** 


u 


over the whole layer. A relation between the distributions of enthalpy and shear 
stress is then found from the energy equation, using a turbulent Prandtl number a 
which is assumed constant across the layer to relate eddy conductivity to eddy 
viscosity. The final expression is similar in form to Crocco’s integral for the 
laminar boundary layer with « taking the place of the laminar Prandtl number co, 
but contains two extra terms proportional respectively to («—o)c, and (a — o)c}, 
which represent the effect of the inner viscous regions. 

The enthalpy integral is evaluated using the stated velocity-shear relation, 
and an expression which agrees well with the available experimental data is found 
for the heat-transfer coefficient as a function of recovery factor and skin-friction 
coefficient. It is also found that the usual quadratic enthalpy-velocity relation, 
exact for a = 0 = 1, remains an acceptable approximation for Prandtl] numbers 
considerably different from unity. 


+ At Cornell University, Ithaca, New York, 1959-60. 
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i. introduction 1.1. Howarth’s transformation 


A number of recent authors on the compressible turbulent boundary layer have 
noticed that, as in the laminar case, the analysis can be greatly simplified by 
using instead of the physical co-ordinate y the variable 


ry 


n=] (p/p) dy (1) 


/0 

introduced by Howarth. Here p is the density at height y above the surface, and 
p, an arbitrarily chosen reference density. For the laminar boundary layer on a 
flat plate, if the product py of density and viscosity is treated as constant, the 
effect of the transformation is to free the velocity profiles of explicit dependence 
on Mach number and Prandtl number. They can be written 

Ulu. = f(9/A), (2) 
where f is a universal function, w,, is the free-stream velocity at the outer edge 
of the boundary layer, and A is the transformed value of a suitably-defined 
boundary-layer thickness 6. 

It seems natural to ask whether any property like this may also hold for 
velocity profiles in the turbulent case. In an attempt to answer this question, a 
set of profiles measured by Lobb et al. (1955) is examined in the present paper. 
These were obtained at Mach numbers between 5 and 8, and momentum-thickness 

Reynolds numbers R, between 5000 and 12,500, with varying amounts of heat 
transfer. When plotted in the form of u/w,, against y/é as in figure 6 the profiles 
are widely different, but the transformation (1) almost eliminates the scatter, so 
that in figure 7 all lie close to a single curve. Except near the inner and outer 
limits, the power law uu, = (n/A)"n (3) 


with n = 9 provides a good interpolation to the curve, and from experience of 
incompressible boundary layers one would expect the index 1/n to decrease 
slowly with increasing R, over a wider range. 

This particular form of interpolation has of course no special physical signi- 
ficance, but it does make the calculation of the remaining properties of the mean 
flow fairly straightforward. Most recent workers on the incompressible boundary 
layer have preferred to express the outer velocity distribution in terms of the 
‘defect’ or ‘wake’ laws elegantly formulated by Clauser (1956) and Coles (1956). 
These laws offer more insight into the physical processes at work in trans- 
ferring turbulent energy through successive eddy scales within the boundary layer, 
but in the zero-pressure gradient case they have little advantage in accuracy 
over the simpler type of law used here. Dimensional analysis shows that the 
incompressible boundary-layer profile on a flat plate can depend on two para- 
meters only; in the outer part one of these must be y/d, and no generality is lost 
in using for the second the index 1/n of the interpolating power law, which clearly 
depends on Ry. 

1.2. Law of the wall 
The form of equation (2) is, however, inappropriate in the sublayer and inner 
turbulent core, for we know that in incompressible flow the velocity distribution 
24 Fluid Mech. 8 
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in these regions is given by the ‘law of the wall’. That is, the velocity ratio u/w,, 
where w, is the friction velocity (7,,/p)*, is a universal function of yu,/v. This 
could only be strictly true if the boundary-layer thickness were infinite, but its 
approximate validity at finite Reynolds numbers is established beyond doubt, 
and one is naturally led to try to generalize the law for compressible flows. 

A way of doing this is suggested by the form of the sublayer profile. If the 
friction velocity is defined for a compressible flow in terms of the reference density 
Y To * Potts, (4) 
then, since the shear stress 4(Cu/Cy) in the sublayer is nearly constant and equal to 
T,,, the profile may be written 


= — a YU,/Vo- (5) 


Here v) = (p/P), and it has been assumed as before that py is constant. The 
accuracy of this assumption is discussed in §7. Equation (5) is of the form 

3 = g(yu,/Vo) (6) 
and by analogy with the incompressible case one might expect an expression of 
this form to hold also in the buffer layer and in the turbulent core. To agree with 
the well-known expression for incompressible flow, g would have to approach the 
form A log (yu,/v))+ B in the fully turbulent core, i.e. for yu,/v, > 30 say. If, in 
addition, the physically realistic requirement of a continuous derivative, i.e. of 
continuous viscous shear stress, at the intersection with the sublayer profile (5) 
is imposed, we are led to the compressible version of an expression derived by 
Squire (1948), namely e 


wil log (7 c)+B (7) 


u, Vo 


in which A and B have the accepted values 2-5, 5-5, respectively, and ¢ = 5:3. 
With these values the intersection of (5) and (7) is at yu,/v, = 7:8. 


1.3. Choice of the reference density 


The outer velocity distribution (2) is independent of ,, but that defined by (6) 
does depend on the particular reference density which is chosen. A way of de- 
ciding on the appropriate value is suggested by the ‘intermediate enthalpy’ 
method of correlating compressible and incompressible values of the skin fraction. 
It has been found experimentally—for example, by Matting, Chapman, Nyholm 
& Thomas (1959)—that the wall shearing stress 7,, a distance downstream of 
the virtual origin of turbulence in a compressible boundary layer can be written 
kD Un C/(U,,2/V,,), Where c, is the same function of the Reynolds number as in 
incompressible flow. The suffix m now denotes evaluation of the physical quanti- 
ties at a certain intermediate value h,, of the enthalpy, for which Eckert (1955) has 
given the empirical equation 


h,, = 0-5(h,, +h.) + 0°22(h,—h,), (8) 
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where suffixes w, co and r refer to wall, free stream and recovery (i.e. zero heat 
transfer) values. 

But the velocity distributions (3) and (7) alone contain enough information 
for the skin friction to be deduced, and it is shown in §2 that if the two laws join 
with ou/cey continuous, then 


TylPyt2s = Opie tac |flg) 2? = Di tag a7/¥y) 2, (9) 


where C(n) and D(n) have almost exactly the values used in the recognized 
incompressible power-law formulae. The fact that the power-law formulae for 
skin friction can be derived a priori in this way does not seem to be well known: 
most authors have regarded the formulae, at least for the casesn = 9andn = 11, 
merely as interpolations to more complicated expressions of the Karman-Schoen- 
herr type, which were originally derived by use of a velocity-defect law. But if 
Pp, and fly are the values which apply at Eckert’s intermediate enthalpy h,,, then 
(9) gives precisely the value of skin friction predicted when his law is used to 
generalize the latter formulae to compressible flow. It follows, therefore, that if 
we set 
Po= Pm Fo= lm (10) 
the inner law (6) is consistent with the known variation of skin friction with Mach 
number. These reference values are therefore used in the remainder of the paper. 


1.4. Distributions of shear stress and enthalpy 

Before comparing the velocity laws (3) and (7) directly with the available experi- 
mental profiles, it is useful to calculate the distributions of shear stress and 
enthalpy across the layer on the assumption that they do in fact hold. The shear 
stress is found in §3, using the equations of motion in conjunction with appropri- 
ate forms of the steam function in inner and outer regions. General expressions 
involving only the functional forms of the velocity profile are found first, but 
on inserting the precise distributions (3) and (7) an approximate expression is 
obtained for 7/7,, in terms of w/u,, alone. This is required in §4, where a turbulent 
analogue of Crocco’s integral relating the enthalpy distribution to those of shear 
stress and velocity is derived from the energy equation. 

The derivation calls for some further knowledge of the turbulent transport 
mechanism, in addition to that implied by the distribution of velocity (from 
which those of shear stress and therefore of eddy viscosity followed). For this 
purpose it is assumed, following van Driest (1955) and Rubesin (1953), that the 
ratio of diffusivities of momentum and heat due to the turbulent motion is con- 
stant across the layer. This ratio defines a turbulent Prandtl number «, different 
in general from the laminar value o which corresponds to molecular transport 
mechanisms. Rubesin and van Driest had recourse to mixing length theories to 
provide the distribution of shear stress, but we can now use the distribution 
obtained directly from the velocity profile. Both this and the distribution of 
enthalpy to which it leads are therefore consistent with the equations of con- 
tinuity, momentum and energy, and with each other. This was not the case with 
those obtained by previous authors. 
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1.5. Recovery factor and heat-transfer coefficient 


From the enthalpy distribution, expressions are found in §5 for the recovery 
factor (i.e. the fraction of the kinetic energy of the free stream which is recovered 
as enthalpy at an adiabatic wall) and for the heat-transfer coefficient. In each, 
the more important term depends only on «, but correction terms proportional 
to (4-0) (p.,¢,/p,,) and to (2—0@)(p.¢,/p,,)*, respectively, appear. The value of 
a cannot be obtained directly from experiments, but it is possible to infer a value 
from that of the recovery factor, which can be measured. The value of « so found 
can then be used to predict the heat-transfer coefficient which corresponds toa 
particular recovery factor. The resulting heat-transfer coefficients lie well within 
the range which would have been expected from experiments. 

Finally, in §6, the experimental velocity profiles referred to are examined to 
show their compatibility with the model of the boundary layer which has been 
elaborated in the paper. 


2. Expression for skin friction 

As remarked in §1.3, an expression for skin friction can be derived directly 
from those for the inner and outer velocity profiles, on the assumption that the 
gradient Cu/ey is continuous at their intersection. Using (2), (6) and writing 
nu,/Vo = €, we should have at this point 


© Site) Mh oy 
U,, =A(7s] ses g(c), 


L ou Ve hv re 
7a 

U,06 u,A° \u.A Ux 

Eliminating ¢ between these two equations clearly yields a relation between 

u,/u, and u,A/v, or, what is the same thing, between (w,/u.,)? and u,A/v). In 

particular, the power law (3) can be rewritten 


uN u —(n+)Din /y A ln > c 
= | "| | : Cin = Kin (11) 
U. Us Vo 
say, whence 
2 2(n+1) 
7 uU, \* a w-A\ = 
w= (—t) = K-2ntntn (=) (12) 
Po tx U, Vo 


Applying the second condition that cu/ef is continuous at the intersection of 
(11) with the logarithmic law (7) (from which ¢ can be omitted since the join 
occurs at values of € much greater than 5-3), A is found explicitly as 


K = nAexp[(B/nA)—1]. (13) 


The values of A for n = 7, 9, 11 are tabulated in table 1, and the corresponding 
outer curves are plotted against Cin figure 1. The fact that the power and logarith- 
mic laws in each case lie close together for a fair distance on either side of the 


point of contact shows the joining condition to have been a realistic one. 








fo 


2 


sl] 





n 


f 
n 


— 


Velocity distribution in the compressible turbulent boundary layer 373 


The momentum thickness may be evaluated using the power law (3), since 
integration smooths out errors produced by the departure of the profile from this 
form at the inner and outer edges; we obtain 





3 
O=]| (pu/p.u,)(1—u/u,.)dy = {n/(n+ 1) (n+ 2)} (pp A/p.). (14) 
/0 
n K 2/(n +1) C Range of Ry 
7 8-82 1 0-0124 0-:0130 100-5000 
9 10-57 3 0-0087 0-0088 500—50,000 
11 12-36 1 0-00645 0-00655 3000-600,000 
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FicurE 1. Various ‘inner’ velocity laws. 


Xeplacing A in (12) by this value of 6 yields the first of the expressions (9) for 
skin friction. The second, in terms of a Reynolds number based on z, follows 
it by writing the integral momentum equation as 

A( pot! n)/d(Unt]V,) = (u,/U.)” 

and setting 6 = 0 when x = 0 after integration. The constants C and D are 

related to K by 
(n+1)O8"+D = (n +3) D¥"+9 = n(n+ 2)1K-. 

The values of C obtained using K from (13) in this relation are listed in table 1, 
together with values C’ corresponding to the accepted incompressible power-law 
expressions 3c, = C’R,*"*) given in the three cases by von Karman (1921), 
Young (1953) and Falkner (1943). The ranges of applicability of the incompres- 
sible laws are those suggested by Culick & Hill (1958). 
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The agreement between the value of C found by the method just described 
and the accepted value C’ is close for the 4th and 4th power-laws, and within 5% 
for the quarter power-law which holds at the lowest Reynolds number. With 
slightly different values of K (8-85, 10-50, 12-23, respectively) C could be made to 
agree exactly with C’. The corresponding power-law velocity curves would still 
approximate closely to the logarithmic expression over a considerable range of 
R, and might just as well have been chosen in the first place instead of those given 
by (13). 


3. Shear-stress distribution 
The equations of continuity and motion are 


~ A 


C : 
2 : ra? = 0, 15 
om (Pu) +5, (pv+p'v’) (15) 

Ou —,0u OF P 
pu a (pv+p. ‘by "aS (16) 


respectively, where primes denote turbulent fluctuations, bars the time means of 
their products, and 


ou | we 
T=UM~—puv 
et 


is the shear stress. In view of (15) a stream function yy may be defined by writing 


— 


oy ey — ‘ 
Pm oy = pu, —Pm Cx = pe +p v. (17) 


Different forms for y are appropriate in the inner and outer regions, which we 
now treat in turn. 


3.1. Inner region 


The stream function may be written 


Y= G(T) = 79 G0) (18) 
nn ast ; 
say. Then g(C) = te G'(C) 


and (treating v,, as independent of 2, i.e. for a constant temperature wall) the 
momentum equation (16) becomes 


A 


9 OT 
—pu*A(x) = —, (19) 
cy 
where A(x) = —(d/dx) log (u,/u,.). This is the compressible flow version of an 


expression found in an earlier paper (Spence 1951). On integration 


re 
T| Tw = 1 Pes (Vin/ u,) A(x) [9(f,) ? dg. (20) 
0 


. 


Thus with a knowledge of the universal function g in equation (6) and of the 





Sa 


Ww 
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variation of 7,, with x, the shear stress in the inner part of the layer may be found 
from (19). The result of integrating the logarithmic expression (7) can be written 


a 


t 
| g?dg = fg>{1— Ag+ O(g-*)} 
0 


so that for large enough values of g = u/u,, 


U 


Ti 


Tle = 1= (m/s) Aa“) = 1-Aeyn(*).. (21) 


3.2. Outer region 


The velocity law (2) may be written in terms of a stream function 


yy = u,AF(n/A) = uAF(E) (22) 
say. In this case F(é) = — = F’(&) 


and the momentum equation can be written 


—AF(£) F"(&) = (Pn IPT) 5, - ag (=). (23) 
where A = (p,,u?,/7,,) (dA/dx) may be supposed known. F contains a constant of 
integration which should strictly be evaluated by equating this value for 07/0y to 
that given by (20) at the point of contact of the corresponding velocity laws, but 
as in determining 6/A, an adequate approximation is found simply by applying 
the outer law across the whole layer. Then by (22) F(0) = 0 since y% must be 
constant on £ = 0, and using the power-law velocity distribution we obtain 
F(&) = {n/(n+1)} £11. Therefore, integrating (23), 


1-7 = giitin = (=). (24) 


w Ux) 


The coefficient of £!+?’" here has been set equal to unity in order to make 7/7,, 
vanish at &=1. The coefficient is actually nA/(n+1)(n+2), which equals 
(p.,u2,/7,,) (dO/dx) by (14), so this is a statement of the integral momentum 


equation. 
The expression (24) is not very different from that which would have been 


found by using the inner law in the region where it applies, for logarithmic 
differentiation of (12) shows that 


A(x) = {1/(n + 1) A} (dA/dz), 
whence (21) may be written 


? A yn ([u ) ( *) ( 2 
ele ai —{] 2 * 
: * ae Ta Va. 


w 


This is an approximation to (24) ifn > 1. 
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Finally, (24) may be written in terms of velocity alone as 
? i.) ‘ 
- =1-(—} ; (25) 


Tw Un 


in which form it will be used in the next section to evaluate the enthalpy 
distribution 


4. Enthalpy distribution in terms of velocity and shear stress 


The energy equation may be written 


oh 5,0h og ou 4 
pu—t(pyt+p'v')— =~—4+7T , (26) 
Ox Cy cy oy 
where h is the static enthalpy, and 
k oh _ 
q=7>-—puh 
C, cy 


is the heat flux per unit area towards the wall. In forming equation (26) triple 
products such as u’*v’ have been excluded as is customary. Regarding h and q 
as functions of u, so that ch/Ca = (Cu/ex) (dh/du), ete., and using the equation 
of motion (16), the left-hand side of (26) is seen to equal (é7/éy) (dh/du). After 
dividing through by 7(¢u/cy) the equation may therefore be written 


1dr dh yy 1 dq 
tdudu  tdu 


An eddy viscosity ¢ and an eddy conductivity k may now be defined by means of 


Cu kK Ch 
—pu'v’ =e, —ph'r' =~ —. (28) 
Cc 
Cy p cy 


The laminar and turbulent Prandtl numbers are then 
Y /. Bex is 
c= "uC,|/k, a=eC,/k, 


respectively. In terms of these the ratio of the molecular contributions to q and 


T 18, say, - k ch ou om 1 dh 
Grit, = (c sy) (/ ay) ~ odu’ 


and similarly that of the turbulent contributions (28) is 


ldh 


OnlTm = = 
{ri71 adu 


Combining the last two equations, 


q=41+47 = (72472 dh _ |r (- ! Be 


o a)du \a \o « du 
and the energy equation (27) becomes 
1d {{r _ \dh| ldrdh 
: { -_ a - ==_— .. 29 
saaltet (s ) "1 ial tdudu (29) 
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So far it has not been necessary to make any physical assumptions, and equa- 
tion (29) is exact. To proceed further we may suppose both o and « are constant 
across the layer. The accepted value of ¢ is about 0-7 for air up to about 1000 °K, 
and a value of about 0-85 for a will be inferred in §5 from the known values of 
recovery factor. 


4.1. Integration of the energy equation 
The fact that the viscous contribution to shear stress falls off very rapidly with 
distance from the wall—in fact exponentially with velocity, decreasing almost 
to zero in a region of the turbulent core where the total stress is still not much 
less than 7,,—makes it possible to integrate the last equation in a simple manner. 
After multiplication through by 7!~*, the equation can be written 


d dh d {dh 
(r- “: ) tarts = ee (70, )> (30) 


du du du 


where / = (a/a7)—1. In the region where r; + 0, r may be treated as constant 
and equal to 7,, on the right-hand side. Then integrating 


(i) sof (Eo) 


0 \7., “du du), 
and re-arranging, 


dh ([T = 7 \P-* fe, | 
en ee “ae A(z 1) 2). (32 
du (2) a a] -T— (“) “7 (=) (# oo ) (32) 
The result of integrating (32) with respect to wu, treating 7 again as equal to 7,,, in 
the term involving 7,, may be written, excluding O(f?), as 
‘dh’ , U, yh Fe) 

Me (T) {U)—(e—0) 2 Lig)| ue [K(f) -(a—0) (2) i}, 
where for convenience f and g have been written for w/u,, and u/u,, respectively, 
and 


Naa) (Fy a | 


i eo (34) 
K(f) = | 1 | if, | 
’ ) "7 0 (;) dfs | 0 oa dfs 
g °g 
and L(g) = [ “hdg, N(g) = | "Lg dg. (35) 
0 'w J0'w 


On setting « = o, we are left in (33) with the well-known Crocco integral for the 
enthalpy distribution in the laminar boundary layer at zero pressure gradient 
(which could also have been obtained at once from (32)). To evaluate the last 
two integrals we may set 7, =7,, in the sublayer, and obtain the value of 
T,, = (Cu/oy) in the buffer layer by differentiating (7) with py constant; this 


yields _¢ 
TT, = 1, exp(-75"), 
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according as u/u, = g 5 ¢,( = 7:8). Then for g > ¢,, 


L=L,—£A4, N=NH,- 1 Ag+), (36) 


w w 


where [L,=&,+A = 103, Ny = $1§5+4A? = 56- 


Except in the innermost part of the turbulent core we may set 7,/7,, = 0 and 
simply use the latter values. Values fairly close to these can also be obtained from 
other recognized interpolations for the buffer-layer profile; for instance, Rannie 
(1956) has suggested the formula 

u sl 
a (1/\«,) tanh (,/«;) 
with 1/,/«, = 14-53, for incompressible flow in the region 0 < €< 27-8. The cor- 


responding value of 7,/7,, is 1—K,(u/u,)?, pou by integrating up to the point 


ss 
where this vanishes we find L, = 9-7, N, = 52-8. (The less accurate assumption 
of an —— ne interface w ithout a buffer layer at the point 
where € = A log + B, i.e. at € = 11-7, gives L, = 11-7, Ny = 68-5.) 


4.2. Recovery factor and heat-transfer coefficient 
Since it gives the enthalpy at the outer edge of the boundary-layer equation (32) 
provides expressions for the recovery factor, i.e. the fraction of the kinetic energy 
of the external stream which is recovered as enthalpy at an insulated wall—and 
for the heat-transfer coefficient. With wu = w,, and introducing the skin-friction 
coefficient c, = p,,u7/4p..U%,, (32) becomes 


a 


du 


; ot 
(K(1)—(a—0) N,22%1\ 
"Pm2} 

(37) 
In the case of zero heat transfer (dh/dw),, = 0, and h,, equals the recovery en- 
thalpy h,. The recovery factor is therefore 


= (h,—h,,)/4u2, = 2K(1)- 


(1) —(2-0 7 aes s 1)"| ut 


he—hy =| 4 


Dan = 
Po (4 —o) Noe. (38) 


m 


The inward heat flux at the wall per unit area when h,, + h, is 
- k = = (“) be (R).- ie) 
‘ Cy Cy) w To) » \C 2] du du} 
Elimination of h,, from (37) and (38) shows that the first term on the right-hand 


side of (37) is equal to h,—h,,. This value may be used to calculate the Stanton 


number, defined by 


1/St = poUalh,—Ny)/Qe = (2/¢,) L(1) — (2px Py)? (a—@) Lg. (39) 
For « = 1 we should have /(1) = 1, and the — hand side of (39) reduces for 
incompressible flow to the form 2/c,—(2/c,)!¢(7) found by Rannie (1956). If 


also a = 1 the equation becomes il St = 4c,, which expresses the Reynolds 
analogy between heat and momentum transfer. In the general case, a Reynolds 
analogy factor may be defined by the ratio St/}c, as calculated from (39). 
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5. Evaluation of the enthalpy distribution 
The distribution of shear stress in terms of velocity as given by (25) may now be 
used to evaluate the functions J(f), K(f) which occur in the enthalpy-velocity 
relationship (33). Writing f for w/w,., (25) is 
T/T, = 1-f?, 


a 






<a & i ] “iE 
2K (f | wae | 
broken lines represent 2/* K(1 broken lines represent a f } 
| | | / 
| _>Values of @ 
- 1-4 
08 | 
0-4 
) | | 
. 0-4 10 











f 


Figure 2. I(f) and 2K(f) for several values of « (p = 11). 


where p = n+2, and J and K, from (34), are 


ry 
I(f)=a| (1-fPytdf, 


/0 


; of hy 
Ki=a[ -syds, [0 -p2dh, 


0 

These functions have been computed for p = 9, 11 and a = 0-5(0-1)1-5. Z(f)is an 
incomplete beta function, but the standard tables do not cover the required 
range of « and p; the functions were therefore found by quadrature after writing 


fs = 1—a 
(fps tmz, (FR) =i), (AN -Pyy = ke). 


(Pah) : : 1 raf) : Ps, 
i(z)dz,, K(f)= | 0(z,) dz, | k(z.) dz», 
/~0 


Then for « < 1 


I(f) = 


and for a > 1 





~0 0 


lid an 


af) raf) 2 
I(f) = 5 : | k(z,)dz,, K(f) = (2 = | k(e,) da | (22) dz. 
2 F 2 


) 0 0 
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A check is provided by the fact that 
i ] 
-_ a — 
10) =a! ( (2+, if. (40) 


which can be found from tables. 


“| 











FiGuRE 3. The values of /(1) for p = 9, 11, 13, and 2K(1)/I(1) for p = 11. 


Figure 2 shows the results for x = 0-6-1-4 and p = 11 (i.e. for the 3th power 
velocity profile). The broken lines representing af and f2K(1) lie close in all cases 
to the curves of /(f) and K(f), respectively, and it follows that within the limits 
imposed by ignoring the terms involving c, in (33), the enthalpy variation across 
the boundary layer can be well approximated by a quadratic in velocity, i.e. by 


u i 
h(u) = hy + (hy lap) = — trv? (41) 
x 
over the range of Prandtl numbers (or recovery factors) covered by the computa- 
tion. 





Ee 





















0) 





Figure 3 shows the variation with « of /(1). The curves for p = 9 and p = 11 lie 
very close together. A curve of 2K(1)/I(1) is also shown; in this case the results 
for p = 9 and p = 11 were indistinguishable. The variation of 2K(1)/Z(1) is fairly 
slow, and its value can be read sufficiently accurately at non-tabular values of « 
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to permit accurate evaluation of K(1), using the exact expression (40) for J(1) 
at such values. A fairly wide range of Reynolds numbers is probably covered by 
the case p = 11; the values of J(1) and K(1) for this case have therefore been used 
to express the recovery and Reynolds analogy factors as functions of « for a 
number of values of (w,/w,,)? (i.e. effectively of boundary layer Reynolds number 


Pox 4/u,,) in figure 4. 
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As was stated in §1.5, values of « cannot be found directly from experiment. 
We do know, however, that for air at Mach numbers up to 8 the recovery factor r 
lies somewhere between 0-88 and 0-91. The curves of figure 4 have therefore been 
cross-plotted in figure 5 for constant values of r in this range, to show the variation 
of the Reynolds analogy factor with (w,/w,,)?. The actual magnitude of the factor 
St/}c, and its decrease when c, decreases, i.e. when pu. 9/,, increases, are in 
good agreement with those found experimentally, which have been reviewed by 
Seiff (1954). The values of « in the cross-plot are in the range 0-85-0-9. 


6. Experimental velocity profiles 
6.1. Outer law 

Figure 6 shows a set of profiles of u/u,, against y/d measured by Lobb et al. (1955) 
at Mach numbers between 4-93 and 8-18 and stagnation temperatures near 
300 °K. The M = 8 profiles lie near a }th power law, and those for M = 5neara jth 
power law. The effect of the Howarth transformation is illustrated in figure 7 by 
plotting the profiles with 7/A as the ordinate. This was found in terms of y/@ by 
writing 


2 


” 
y=) (hfh,,)dy (42) 
~ 0 
and using the quadratic enthalpy-velocity relationship (41) with w/u,, = (y/A)!". 
Integrating and inserting @ from (14), the relation 


n/A = (y/O) [a+ b(u/u,) —e(u/u,)*]- (43) 
is obtained, where 


na = (n+1)(n+2)(h,,/h,), (n+2)b = (h,—-h,)/h,, (n+ 1)e = (h,/h,)-1. 


The experimental values of y/, u/u,, and a, b, c (with n = 9) were used in (43). 
As plotted in figure 7 the profiles have collapsed almost entirely to a single curve 
lying close to the }th power law. This justifies the use that was made of the law 
in simplifying the integration with regard to 7. 

To avoid overcrowding the figure, only five of the profiles measured by Lobb 
et al. are shown. The same collapse to a single curve was, however, obtained with 
the remaining nine profiles as well. 


6.2. Law of the wall 
The inner parts of the same set of velocity profiles are plotted in the form of 
u/u, against yu,/v,, in figure 8, which also contains data from tests at the R.A.E. 
by Monaghan & Cooke (1953). The calculation of 7 from y by means of (42) is a 
little more complicated in this case, since it is necessary to take account of the 
dependence of the enthalpy distribution in the sublayer and buffer layer on 
T,/T,,, a8 represented in (33) and (36). The details need not be given here. 

The friction velocity was found from the experimental values of p,,u,,4/1,,, 
using the }th power-law formula for skin fraction in the form given by (9). The 
intermediate enthalpy was found from the wall and free-stream values, using 
Eckert’s formula (8) and assuming a recovery factor r = 0-89. 
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Ficure 6. N.O.L. velocity profiles in physical co-ordinates. Key: 
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The values of (w,/u,.)? so found agreed well with those published for the R.A.E. 
data, but with the N.O.L. data there were systematic discrepancies of up to 25 % 
in the cases with heat transfer. The values of c,; measured by the N.O.L. decrease 
with wall cooling, which is the opposite tendency to that predicted by the inter- 
mediate enthalpy theory. The measured values were obtained (a) from the 
measured heat-transfer rates, assuming the Reynolds analogy factor to be 1-2, 
and (b) from the slopes of the velocity profiles in the sublayer. Neither method is 
very reliable, and the theoretical values of skin friction corresponding to the 
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9, [Vim 
FicurE 8. Inner behaviour of experimental velocity profiles. ———, 2-5 log (yu,/v,,) + 5:5. 
R.A.E. data 13-4 station: @, with heat transfer; A, zero heat transfer. Other symbols as 
for figures 6 and 7. 


observed local Reynolds numbers were preferred in the present analysis to those 
quoted by Lobb et al., since the bulk of other experimental work tends to support 
the intermediate enthalpy formula. The collapse to a single profile shown in 
figure 8 would not have been as complete if u, had been based on the quoted values 
of c,, although the curves would still lie closer together than when plotted with 
y, V,, and a wall-friction velocity u,, = (7,,/Pr)? in place of 7, v,, and u,, asin their 
paper. 

There is a small but systematic discrepancy between the data from the two 
sources as plotted in figure 8, the N.O.L. points all lying a little above the standard 
curve u/u, = 2-5log (nu,/v,,) + 5-5 and the R.A.E. points somewhat below. This is 
almost certainly due in part to residual Mach number effects not removed by the 


transformation, the presence of which would be likely to show since Monaghan 
25 Fluid Mech. 8 
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& Cooke worked at about M = 2 and Lobb et al. at hypersonic speeds. The dis- 
crepancy may also reflect the fact that even to-day there is no general agreement 
among experimental workers on the precise values of the constants A, B in the 
logarithmic law (see, for example, Table 1.1 of the forthcoming volume Incom- 
pressible Aerodynamics edited by Thwaites, in which a number of sets of values 
from the literature are listed). The values A = 2-5, B = 5-5 used in this paper are 
probably the best known. The reason for this measure of disagreement may lie 
in the fact that some dependence on a Reynolds number such as w,, 6/v,, should 
strictly be included in the inner law; but in part no doubt it also arises from the 
difficulty of measuring skin friction accurately. 

In a recent note Yasuhara (1959) has attempted a correlation of some experi- 


mental profiles in terms of a defect law of the form (u,,—w)/u, as a function of 


n/A. The fact that this was not very successful may also have been due to un- 
certainty about the skin friction. 


7. Concluding remarks 


From the evidence of figures 6 and 7 there seems little doubt that the method 
of collapsing the outer four-fifths or so of the velocity profiles as a function of 7/A 
is valid, and this justifies a posteriori the calculations of the shear stress and 
enthalpy distributions in §$3 and 4. Other investigators have remarked on the 
fact that the index n in the approximating power law (y/d)!” decreases with 
Mach number—as shown in figure 6—but this change in n is conveniently re- 
moved by the present transformation. This result, and the further fact that the 
enthalpy-velocity relation is approximately quadratic, can be used empirically 
to treat more complicated boundary-layer flows as well: for instance, the relation 
between the compressible and incompressible values H and H,, say, of the form 
parameter 6*/4, is found from (2) and (41) to be 


H = (h,,/ho)H,+(h,/h.)—-1, (44) 


and this equation has been used (Spence 1960) to calculate the growth of the 
boundary layer in a pressure gradient by an extension of the methods available 
for incompressible flow. 

The uncertainty about the skin friction in the N.O.L. tests make the validity 
of the inner law (7) less certain than that of the outer law. One can certainly draw 
the more limited conclusion that if uw, is defined from: the value of 7,,, predicted by 
equation (9), irrespective of whether this is correct, then w/w, is somewhere near 
the standard logarithmic function of yu,/v,, given by (7). Thison its own, however, 
is not a very tidy result, for if p,,u? so defined were not in fact equal to the wall 
shearing stress, the profile in the sublayer could not be written wu/u, = yu,/v,, asin 
equation (5), and the reasoning leading to the functional form u/u, = g(nu,/V,,) 
in the logarithmic region would no longer apply. 

It was also assumed in deriving this expression that the density-viscosity law 
could be approximated for constant pressure by pu = constant, i.e. excluding 
effects of dissociation by oc T. This might at first sight seem to rule out hyper- 


sonic flows at the extreme stagnation enthalpies typical of re-entry. But the 
relation can always be taken in this form over a suitably restricted range of 
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temperatures, and as table 2 shows, the width of the range over which a single 
constant of proportionality could be used becomes surprisingly wide at high 
temperatures. The table is constructed from the results for dry air calculated by 
Hirschfelder e¢ al. (1955). The assumption is clearly satisfactory for the tem- 
peratures between 100 and 300°K covered by both sets of data analysed, and 


T (°K) 100 200 300 500 1000 2000 3000 5000 
h/T x 10° 7-02 6-68 6-17 5:36 4-26 3°34 2-90 2-41 


(z/em see °K) 
TABLE 2 


should work moderately well, at least across the sublayer, to which the viscous 

effects are mainly confined even for wall temperature as high as 1000 °K. Since 

the viscosity enters the expression for skin friction only to a small power, this 
will not be greatly in error. But the use of the intermediate enthalpy formula 
would of course be highly speculative at such temperatures, and departures from 
the inner law might well become more marked. On the other hand, there is no 
particular reason to expect the outer law to break down at higher temperatures, 
since it does not depend on viscosity; but the conclusions drawn from it regarding 
the enthalpy distribution might be slightly modified by the variation of specific 
heat and Prandtl number. 
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Flow between torsionally oscillating disks 


By S. ROSENBLAT 


Department of Mathematics, King’s College, London 
(Received 4 November 1959) 


Two parallel infinite plane disks, between which is contained a viscous fluid, 
oscillate torsionally about a common axis. Of specific interest are the cases 
where (i) one disk only is in motion, and (ii) the disks oscillate 180° out of phase, 
but with the same frequency and amplitude. The basic parameter is found to be 
the Reynolds number derived from the frequency, the kinematic viscosity and 
the disk separation. Solutions of the Navier-Stokes equations for both small and 
large Reynolds numbers are developed and the transverse and radial-axial flows 
investigated. 


1. Introduction 

In a recent paper (Rosenblat 1959) the author examined the flows resulting 
from the small torsional oscillations of an infinite disk in a viscous fluid otherwise 
unbounded and at rest. It was found that linearizing expansions depending on 
the smallness of the oscillation amplitude were insufficient to describe one of the 
velocity components, namely, the mean steady radial-axial component. This was 
essentially because the method does not allow for the fact that non-linear inertia 
terms predominate outside a thin boundary layer, though they are negligible 
within it. 

As a sequel, the present paper investigates the analogous problem in which the 
fluid is now bounded by a second, parallel disk at a given distance. Two cases only 
are studied: (i) when one disk performs small torsional oscillations and the other 
is at rest; and (ii) when both disks oscillate with the same amplitude and fre- 
quency, though with a phase difference of 180°. 

Commencing with the Navier-Stokes equations, we seek solutions on the 
assumption that certain non-linear inertia terms may be omitted. It is found that, 
unlike in the above-mentioned problem, such solutions are in fact obtainable 
with all boundary conditions satisfied. Moreover, it is subsequently shown that 
the approximate method is valid provided a single condition relating the ampli- 
tude of the oscillation with the Reynolds number of the flow holds. 

Although exact solutions can be developed throughout they prove cumbersome, 
and approximations for low and high Reynolds number are more convenient. The 
former apply for R < 10incase(i)and R < 40 incase (ii); thelatter for R > 20 and 
R > 80, respectively. 

Transverse velocities are first calculated and are seen to have a boundary-layer 
character at large R. The rotational motion of the disks gives rise to centrifugal 
forces which in turn cause the fluid to be thrown radially outwards. Hence 
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a radial-axial flow is set up which has a mean steady component as well as a fluc- 
tuating component. The behaviour of the steady term is of particular interest and 
is fully discussed. 

In order to balance the fluid expelled by centrifugal action a radial inflow 
between the disks is required. Hence a radial pressure gradient is induced which 
maintains this flow. It is found from the equations of motion to act inwards and 
to be constant. At high Reynolds numbers it is the dominant factor outside the 
boundary layers so that the flow there has a Poiseuille-type profile. 

The case of large amplitude oscillations is not considered here. In the limit as 
the amplitude becomes infinite the problem becomes that studied by Batchelor 
(1951) and Stewartson (1953), where two parallel disks rotate steadily about 
a common axis. Some formal comparison of results may therefore be made. 


2. One disk oscillating 

Consider a body of fluid bounded by two parallel disks which are represented 
by the planes z = 0 and z = d in a cylindrical polar co-ordinate system. The axis 
r = 0, perpendicular to the disks, is the axis about which they may perform 
torsional oscillations. If vu, vand w be respectively the radial, transverse and axial 
velocity components, p the pressure, p the density and v the kinematic viscosity, 
then the Navier-Stokes equations of motion are 


A A A 9 ~ AD ny a9 
ou ou ou vs 1 Cp Cu lew wu cu 
or a 00 go — 8 — ee on es (1) 
ct cr ca fr per Ore ror fr © 
cv =v cv ww Oy ld vi dw - 
a 6-10 — = Pl be al, (2 
ct cr ca rT OC LCF £* Cz 
Cw 00w Ow 1 Op Cw low cw ; 
—+uU—+w- =—--~— +155 <- tals (3) 
ct cr Oz pc CF FOF Cz) 
while the equation of continuity is 
0 ow 
~(ru)+— =0 (4) 
er 0z 


We take the disk z = 0 to perform torsional oscillations of frequency and 
angular speed (2, while the disk z = d remains at rest. The boundary conditions 


are then 


uw=w=0, v=rQe'™ on z= 0,) 


| (5) 


> 


e=v=w=0 on z=d. 
Equations (1)—(5) can be satisfied by writing 


rQ2eF 2 
ae » = —2d ~~ Fly,7), 
-— y,T), U d- (y,7),| 


ple = Q2d?Ply,7)+402r?K(7), y=zid, r=nt. 


v= rQe'g(y), wu (6) 


That is, we seek a solution in which the radial and transverse velocities are linear 
functions of r. The equations of motion then establish the required form of the 
pressure, with P (y,7) and K(r) unknown functions of the variables indicated. 
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On substitution from (6), equations (1) and (2) become 


02k CF \? oF ; 1 BF 
<<. ae = i) ee Sree 7 
dy or -+(; ) {(e) ” | i as > *. “) 
5 Q)\?27 oF dg 1 dg 
9 — A & 
and ig+2 | = E ay ra Ray?’ (8) 
where R = nd?/v (9) 


is the Reynolds number of the flow. The boundary conditions (5) become 


._ 
F= ,=0, g=1 on y=) 
cy 
i ! (10) 
af a 
F=—=g=0 on y=1. | 
cy 
Kquation (3) reduces to 
oF oP 2 eF 
ee - 1] 
5 ee ey Rey? sie 


and serves merely to determine the axial pressure gradient. 


Transverse component 


On the assumption that the amplitude of the oscillations, i.e. Q/n, is small, it 
appears permissible to neglect to first-order the non-linear convection terms in (7), 
(8) and (11), though retaining the non-linear centrifugal term (ge'’)?. (Some 
remarks concerning the validity of this assumption will be made later.) For the 
transverse flow we then have to solve 


d’q 
dy? ,—tkg = 0 


with g(0) = 1 and g(1) = 0. The solution is 


sinh ,/(tR) (1—y) 
= =, 12 
I sinh ,/(72) me) 


Equation (12) leads to, in real notation, 


[cos ,/(4R) y cosh ,/(4R) (2 —y) — cos /(4R) (2 —y) cosh ./(4R) y] cos nt 
vo + [sin /(4R) ysinh /(4R) (2—y) —sin (42) (2—y) sinh / (42) y] sin nt 
rQ cosh , (2) — cos ,/ (22) ; 
(13) 
Expanding this for small values of the Reynolds number, we get 
5 = (l-y) iC _ = 8 + Sy — 12y? + 3y9)| cos nt + ny (2—y)sin nt + O(R?), 
. (14) 
with amplitude approximately 


|r| 


al 
~ (1-y) [ - Fay 2-9) (2-2y +99) (15) 
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and phase angle 


tan! 1 Ry(2—y). (16) 

For large R, (13) becomes 
v . ‘ 
=O ~ e-VGak ” cos (nt —./(4R) y), (17) 


the well-known shear-layer solution, whose magnitude is 


[2 

—! — e-ViRR)y. 18 

rQ - 
Numerical evaluation of |v|/rQ from (13), (15) and (18) for various values of R 

indicates that the approximation (15) is valid for R < 10, and the asymptotic 

form (18) for R > 20. Curves illustrating |v|/rQfor R = 5, 50 and 200 are given in 

figure 1, the first derived from (15), the two latter from (18). 


10 











Jv] /rQ 


Figure 1. One disk oscillating: amplitude of transverse velocity |v|/rQ for R = 5, 50, 200. 


We see that for small Reynolds numbers the fluid has transverse motion for all 
0 <y< 1. Inthelimit R — 0, the velocity profile is linear in y (the broken line in 
figure 1) with zero-phase angle. For non-zero R, the profile assumes the form of 
a polynomial in y, and the fluid acquires a phase lag with respect to the oscilla- 
tions of the plate. The curve R = 5 shows that |v|/rQ decreases with increasing R; 
at the same time this decrease is less for small y than for large y. That is, there is 
already a trend towards concentration of the moving fluid near the oscillating 





plate. 
At about R = 50, a boundary layer adjacent to the rotating disk becomes 


apparent, and its thickness has order of magnitude 


1 1 /p d |v (19) 
y~ = or z~ = ‘ ‘ 
y JR dun JR Nn 
Thus for large & the presence of the stationary disk has no effect on the transverse 


motion whose velocity profile is exactly that in the case of a single disk in an 
otherwise unbounded fluid. 





392 S. Rosenblat 


From equation (14), the skin friction for small R on the disk z = 0 is 


‘a _— R? R : 
fe C a. ‘all 1+ =) cos nt 3 sin nt (20) 
and on the disk z = d is 
Cv a= 7TR2 | R ; % 
n(- “a 7 |(1— ggg) Cosme +5 sin a ' (21) 


When F + 0, the phi forces are equal on the two disks. As R increases, the 
shear on the rotating disk increases in magnitude, while that on the stationary 
disk decreases. 

Similarly for large R, from equation (17) 


Q /R 
. (; A. -- 5 (cos nt — sin nt) (22) 
and Hl - Me a (2R) e-Va® ~ 0, (23) 
Oz z=d d V = 


Steady radial-axial component 
Reverting to equation (7), we see that the radial-axial flow will be given by 
oF 1 @F 
—(ge'’)? = —K(r)+— —, 24 
We = — Kit aa (24) 


provided the non-linear convective terms can again be neglected. From (12) we 


btai 
bails (get)? 1 [cosh Ay, —cosAy,] _ 1 feosh(1+7)Ay,—1 gait (25) 
ge “= = -~ 
J 2} coshA—cosA |~ 2| cosh(1+i)A—1 
where A=J(2R), y,=1-y. 
This suggests that there are solutions of the form 
F(y,7) = fly) +h(y) "7 (26) 
and K(r) = K,+K,e™. (27) 
Substitution of (25)-(27) into (24) yields 
2 d*f ee. cosh Ay, — cos Ay, (28) 
A2 dy ° 2] coshA-—cosa 7 
poe _— i dh a 1 feosh (1 +i)Aay, —] (29) 
A? dys dy 2} cosh(1+7)A—1 | 
The boundary conditions are 
P=f’=0 on y=0 and y=1.,) 
f=f y (30) 


h=h'=0 on y=0 and y= 1.J 


The above equations indicate that the radial-axial flow has a mean steady 
component and a fluctuating component of frequency twice that of the oscillating 
plate. The radial pressure gradient is similarly compounded of two such terms. 
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The solution of (28) for the mean steady component is found to be 


1 ] 1 
= Si + Si ~ sinh A + si — 3y?+2 
fly) Sao B (sinh Ay, + sin Ay,) a (sinh A + sin A) (1 — 3y? + 2y%) 


+ }(cosh A + cos A) y(1 — 2y + y”) — 4y°(1 -»). (31) 


sinh A + sin A) y(1—y) 


, 1 3 
Fly) = (cosh A — cos A) |- sien + eon 2 | 


+ 4(cosh A + cos A) (1 — 4y + 3y”) — dy(2- 3y)| (32) 


(33) 


; : 3 [2A+A(cosh A + cos A) — 2(sinh A + sin A) 
with K,= , 


A3 
With A and y, replaced by ,/(2) and 1 — y, these expressions reduce to, for small 
Reynolds numbers, 


cosh A — cos A 


R 9 9 /« . 
fe (99 -y (B-y) + (RY), (34) 
7 R —_— ae ‘ = 
PF Faq ylt —y) (6— 15y + Sy?) + O(R) (35) 
and Ky © 35 + O(R?). (36) 
For large R the corresponding approximations are 
] | 
Peal = 2 a 2 2 —V 2R)y. a7 
f~ yl pgp (l-9  + 29) + pope (37) 
, ; 3 init ve 
f'~}l-y(- 3y) +5 ; ary wv" —y)de-very (38) 
. £2 3 ' 
and Ko~ a(5-7am): (39) 


These solutions are illustrated in figures 2 and 3. Figure 2 depicts schematically 
a typical streamline of the radial-axial steady flow at each of the selected Rey- 
nolds numbers Rk = 5, 50, 200; figure 3 shows the dimensionless radial velocity 


f’ = nu/rQ?. When R = 5, f’ = 0 at y = 0-47 and the curves for y < 0:47 and 
: y y 


y > 0-47 have almost the same shape. As RF increases, however, there is a marked 
departure from this symmetry. When R = 200 we see that wu is positive for 
y < 0-38, where both centrifugal force and pressure gradient are active. In this 
region the exponential term in (38) is of the same order of magnitude as the first 
term and the profile is of the boundary-layer type. For y > 0-38 it is seen from 
figure 1 that the centrifugal force is vanishingly small, and the flow near the 
stationary disk behaves essentially like Poiseuille flow with the pressure gradient 
(39). 

It is important to note here that a solution for the steady radial-axial flow has 
been obtainable in which inertia terms are neglected, and yet all the boundary 
conditions are satisfied. Such a solution could not be found in the single-disk 
problem (Rosenblat 1959). An examination of the behaviour at high Reynolds 
numbers throws light on the contrasting features of the twe problems. 
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In both cases centrifugal force, acting within a boundary layer of thickness 
J(v/n), causes fluid to be thrown radially outwards. As compensation, fluid is 
drawn inwards along the axis towards the moving disk. In the single-disk case, 
the quantity of fluid required for this purpose can be obtained in a purely axial 
motion from, in effect, a source at infinity whose strength is determined by 
continuity considerations. Hence no radial flow takes place at large z, and there 














0 40 


FicureE 2. One disk oscillating: typical streamlines of steady radial-axial flow 
for R = 5, 50, 200. 
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vi 
FiGuRE 3. One disk oscillating: steady radial velocity f’ = nu,/rQ? 
for R = 5, 50, 200. 


is no radial pressure gradient. When a second disk is present, however, the 
condition w = 0 on z =d has to be satisfied. Consequently the fluid thrown 
radially outwards near the oscillating disk must be balanced by fluid sucked 
radially inwards. To maintain this inward flow, a radial pressure gradient is 
induced. 

Again in the single-disk problem the centrifugal force is zero outside the 
boundary layer. The hitherto-neglected inertia terms become comparable in 
magnitude with the viscous term and hence important. The result is the formation 





h 


h’ 
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of a second boundary layer within which the steady radial flow takes place. In 
the present problem, on the other hand, although outside the boundary layer the 
non-linear inertia terms are again large compared with the vanishingly small 
centrifugal force, they are still negligible in relation to the radial pressure gradient 
which now dominates (in general, see later). Hence the boundary conditions are 
satisfied and there is no second layer formed. 


Fluctuating radial-axial flow 


Equation (29) yields the time-dependent component of the radial-axial flow. For 
low Reynolds numbers the solution is found to be 


R iR 
b = —— y®X(1 —y)2 | (3—y) ——— (29 + 167y — 1504? + 3043 8), 
h = 5559-9) (s Y) — 75 (29 + 167y — 150y? + 30y )] +0128 (40) 
, R i 
ee iR. 7 > wee ] 
x (6 l5y + 5y?) — 190 (°8 + 385y — 143y? + 1050y? — 210y")| +O(R%) (41) 
- 3 109R 
ag tien 2 9 
and K,= 50 ~ 4300 +0O(R?). (42) 
For high R we obtain 
; (./2 — 2) (.(2iR) — 1) (2—,/(2iR)) _,.,, 
8i(2—,/(2iR))h ~ sb ~(J/2+%e+ —> e-2V Ry 
_ (2 y (26k) —2— V2) | veimy 4 ( 2-2) Vveima—» 
+. ; - o—V (et —_—_—__— ev u ~ 4 
J(2iR) ; +" [(iR) * : 
8i(2—.(2iR))h’ ~ 2—,/2 —2(2—,/(2iR)) e-2¥ Py — (2. ,/(2iR) — 2 — (2) e-V AB 
as (, ie 2) e-v (27R)(1—y) (44) 
ae LD. 
and K,~ =o (45) 


4(,/(2iR) — 2) 

The principal feature of this component of the flow is that at high Reynolds 
numbers a boundary layer forms near the stationary disk as well as near the 
rotating disk—in contrast to the steady component. This arises from the last 
term in equation (44). When y is nearly 1, (44) reduces to 

| (y2—2) [l—e-Veina-v] 
8i(,/(2¢R) — 2) ' 
a typical boundary-layer profile, which is very similar to the solution for flow in 
a pipe under a periodic pressure gradient (Schlichting 1955, p. 199). In both cases 
the formation of a boundary layer is due to the fact that the vorticity diffusing 
away from the wall is constantly changing in sign. 

Outside the two layers the velocity—given then by the first term of (44)—is 

independent of y, since only the constant pressure gradient is there active. 


Remarks 
Finally, we examine the permissibility of neglecting the inertia terms of 
equations (7) and (8), with particular reference to flow at high Reynolds numbers. 
From (38) and (39) it is seen that for large R, f’ = O(1) and K, = O(1/R). Therefore 
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the term (Q/n)? f’? which occurs in equation (7) would appear to be small com- 
pared with the linear terms if 
(Q/n)? < 1/R. 

Since 1/,/R = (1/d)./(v/n), and Q/n is the amplitude of the disk’s oscillations, 
we can say that equation (28) for f and f’ is a valid approximation provided that 
boundary-layer thickness 

distance between disks ’ | 


efle| 
n n 


An alternative way of expressing this condition is as follows. Since by equa- 
tion (6) the steady radial component is u, = (rQ?/n) f’, the rate of convection by 
the radial flow through a typical distance ris Q?/n. Also the rate at which vorticity 
diffuses through a distance d is v/d?. Hence (46) is equivalent to 


amplitude of oscillations < 
(46) 


rate of convection by wu, < rate of diffusion through distance d, ) 


(47) 
i.e. OQ2/n < v/d?. 

The right-hand side of (47) represents the influence of the second disk. If this is 
absent, d > «, v/d? > 0 so that convection is always significant—analogous to 
the result obtained in the single-disk problem. In the present case inertia need be 
included only if (46) or (47) is not satisfied. We would then have to solve a ecompli- 
cated non-linear equation very similar to that encountered by Batchelor (1951) 
and Stewartson (1953) for the steadily rotating disks. 

These difficulties arise only for the steady radial-axial component of velocity. 
For it is quickly seen that linearization of the equations determining the trans- 
verse and fluctuating radial components does not require such a stringent condi- 
tion. In these cases it is in fact only necessary that the amplitude Q/n be small 
compared with unity. 


3. Both disks oscillating 
We now consider the case where the two disks are oscillating with the same 
frequency and angular speed but 180° out of phase. This is analogous to the 
problem, studied by Stewartson, of disks rotating steadily with the same speeds 
in opposite directions. As in that problem we anticipate that the results will be 
symmetrical (or anti-symmetrical) about the plane y = }. 
The boundary conditions (5) are replaced by 
u=w=0, v= rQetr on z=0,) 
; (5a) 
u=w=0, v=—-rQe™ on z=d.J 
Kquations (1)-(4) may again be transformed using (6) and for the transverse flow 
lead to dq 


—iRg = 0, 
dy? iRg 


as previously, but with g(0) = 1 and g(1) = — 1. The solution now is 
_ sinh /(¢R)(1—y)—sinh y(¢R) y 
; sinh (7) ‘ 


(12a) 
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For small Reynolds numbers this yields a velocity profile 


: (1 —2y) | ! = y(l—y)(1+3y—- ay) cos nt + yl —y) sin nt +0(R°). 


rQ- ~ 360° 
(14a) 

Hence the amplitude is given by 

lv| 2 

od | ~ _ — 2y + 2y? 5 

a} 2y)|1 360 Yl — 9) (L— 2y + 2y ] (15a) 
and the phase lag is tan! 1 Ry(1—y). (16a 

p 5 6 Yy Py 


The skin friction is given by 

cv Ov 2urQ rR | + 
(=) = 1(; wie (1+ ) cos nt — sinnt|. (20, 21a) 
02] -=0 02) <a d 12 


720 


These approximations are found to be valid for about R < 40. The amplitude 
profile for R = 20, illustrated in figure 4, is antisymmetrical about the plane 


Seon i 











}v|/rQ 


FIGURE 4. Two disks oscillating: amplitude of transverse velocity |v|/rQ 
for R = 20, 200, 800. 


y = 4in which there is zero transverse flow. The shearing forces on the two disks 
are equal in magnitude, though they will be opposite in sign. 
The asymptotic solution, useful for R > 80, is 


v 


a e-VQGR)¥ cos [nt — ,/(4R) y] —e-v #4-W cos [nt — /(4R)(1—y)], (17a) 


and this leads to 


7 (= = (<) ~ _ Ht (BR) (cos nt —sin nt). (22, 23a) 
2) ent 02) <a d 

Figure 4 shows also the velocity amplitude for R = 200, 800 derived from (17a). 
It is seen that a boundary layer of thickness z ~ ,/(v/n) forms adjacent to both of 
the disks. The expression for the skin friction is identical with (22), when one disk 
is stationary. 

For the radial-axial flow equation (24) again holds, with the same boundary 
conditions. Neglect of the inertia terms is allowable here also provided the 
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condition (46) is satisfied. The substitutions (26) and (27) yield two equations 
analogous to (28) and (29), with a difference in the last term in each owing to the 
changed value of (g e'7)?. 

After some tedious algebra we obtain the following results for the mean steady 
component of the flow. For small R, 











— 
f= 0% |! —y)? (1—2y) + O(R?), (34a) 
4; R : ee 23 or 
f= soy(t—y) (1 — Sy + Sy?) + OCR) (35a) 
and Ky = 7p + O(R?). (36a) 
————________ 
| "i 
10 pees. = = iw 
R = 800/2= 99, e Pei. 
[i oh 
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i ~- | 057 
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FIGURE 6. Two disks oscillating: 
steady radial velocity f’ =nu,/rQ? 
for R = 20, 200, 800. 


Ficure 5. Two disks oscillating: typical streamlines 
of steady radial-axial flow for R = 20, 200, 800. 


For large Reynolds numbers, 


l 


. l 9 D 2k A 
f~ ty —9) (1-29) — ppp (1 —2u) (1 + 2y— 28) + Zapp [eV OY — eV O—M, 


j (374) 


3 
y(l—y)—1fe-VeRy e7¥ (2R)(1—y) 38a 
(oR)! y)— Fl ] (38a) 


f’ ~ 4(1 — 6y + 6y?) + 
A 


. ui, 
and Ky ~ 5 (1- (QR 
(2 


). (39a) 
) 

Figure 5 depicts representative streamlines of this flow at R = 20, 200, 800. As 
might be expected there is symmetry about the plane y = } and no axial flow 
across it. In this sense this plane acts as an interface separating the fluid into two 
self-contained portions, in each of which there is axial flow towards the moving 
disk. The pressure gradient at high R given by (39a) is twice that of § 2, as it has to 
counterbalance the flow in two boundary layers. In both cases, however, this 
pressure gradient is quite small at large Reynolds numbers, which corresponds to 


a result obtained by Stewartson. 
In figure 6 we see the radial velocity component at each of the three values of 
the Reynolds number. In all cases f’ is a minimum (i.e. maximum inflow) when 
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y = }. For R = 20 it is zero on planes y = 0-27 and 0-73. Between these and the 
walls there is radial outflow, due to the predominance of centrifugal action. For 
0:27 < y < 0-73 there is an inflow as the pressure gradient is then the main factor. 

At high Reynolds numbers the profile has the character of a boundary layer 
near both disks, for about y < 0-2 and y > 0-8. In the main body of the fluid its 
shape is described by a quadratic in y, which is well known as the behaviour of 
a fluid under constant pressure gradient. 
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An experimental investigation of the effect of ejecting 
a coolant gas at the nose of a bluff body 


By C. H. E. WARREN 


Xoyal Aircraft Establishment, Farnborough, Hants. 
(Received 29 December 1959) 


An experimental investigation has been made of the effect of ejecting nitrogen 
and helium coolant gases at the nose of a bluff body at a nominal Mach number 
of 5-8. The gases were ejected with ‘swirl’, to encourage them to flow tan- 
gentially to the model surface at ejection, and also straight out. Measurements 
were made of the pressure and temperature on the surface of the model at 
incidences of 0, 4, 8°, and for a range of coolant gas flows. From these measure- 
ments the flow patterns and distributions of heat flux were deduced. 

It was found that ejection with swirl did not in fact lead to an easing of the 
heating problem because the high tangential velocity with which the coolant 
was injected into the boundary layer increased the wall shear stress, and hence, 
by the Reynolds analogy, the heat flux, so that it predominated over the reduced 
driving temperature difference associated with the cooled boundary layer. 

With straight-out ejection it was found that the heat alleviation capabilities 
of the ejected coolant were reduced considerably if the flow rate was sufficiently 
large that the bow shock wave was bulged out. However, provided that the 
external flow is not disturbed, straight-out ejection provides an effective way 
of reducing the heat flux. 


1. Introduction 

A body travelling at high speeds through the atmosphere experiences aero- 
dynamic heating. Ways of alleviating the problems so incurred are accordingly 
of interest. One way is to eject a cool gas from the nose of the body. The gas, 
by flowing back over the body, forms a ‘heat buffer’ between the hot air stream 
and the body. This paper described some experiments that were made to in- 
vestigate some of the fluid dynamic aspects of such a scheme, and, particularly, 
to determine the effect on the local rates of heat transfer. 

In an earlier investigation McMahon (1958) studied the effect of ejecting a 
cool gas from the nose of a blunt body. He tried various methods of ejecting 
the gas, and observed the associated flow patterns. In addition, he measured 
the pressure distributions on the body, and obtained some measurements of 
the overall rates of heat transfer over fairly large segments of the body. Initi- 
ally, McMahon ejected a coolant gas straight out at the nose in a direction 
opposed to the main air stream, but the flow studies and pressure distributions 
indicate severe disturbances to the flow field. It was thought, somewhat 
naturally, that ejection in this manner would not be satisfactory as a means of 




















gen 
iber 
fan- 
ants 
l at 
ure- 


the 
lant 
nee, 
iced 


ties 
itly 
the 
vay 














Ejection of a coolant gas at the nose of a bluff body 401 


reducing the rates of heat transfer. Accordingly deflector caps were fitted over 
the ejection orifice in an attempt to direct the coolant gas tangentially to the 
body surface. The investigation reported here started from this premise. That 
is, the aim in ejection should be to make the gas flow tangentialiv to the body 
surface, and thereby ensure that the inner part of the boundary layer is com- 
posed of cool gas. It was decided, however, that instead of using deflector caps, 
the ejected coolant would be encouraged to flow tangentially to the body 
surface by giving it some swirl in the ejection pipe. It was hoped that the centri- 
fugal effect of the swirl would cause the coolant to flow radially outwards at 
ejection, and that the Coanda effect would subsequently cause it to adhere 
to the surface. The scheme worked, except at the largest rates of coolant ejec- 
tion, but it was found that ejection in this manner is not a satisfactory method 
of easing the heating problem. Although the cool layers of coolant gas reduce 
the driving temperature difference, the high tangential velocity with which the 
coolant is injected into the boundary layer increases the wall stress consider- 
ably, and hence, by the Reynolds analogy, the heat flux. The results indicate 
that the effect of the increased shear predominates, at least for the conditions 
considered here. The matter was not, therefore, pursued further. 

In view of these results attention was directed once more to straight-out 
ejection, for McMahon had subsequently found that this method of ejection 
actually produced an effective ‘blanketing’ of the body. Two flow régimes were 
found. Except at the smaller flow rates the coolant gas emerged as a forward- 
directed jet, causing the bow shock wave to bulge out in much the same way 
as with a solid spike (see, for example, Bogdonoff & Vas (1959)). The normal 
flow pattern is replaced by one in which the main air flow separates ahead of 
the model, thereby leaving a roughly conical ‘dead-air’ region in the vicinity 
of the nose. The flow external to the boundary layer, instead of having a stagna- 
tion point at the nose, now has a ‘stagnation circle’ where the incident flow 
meets the body. McMahon found that the pressures on the body within the 
dead-air region are greatly reduced compared with the pressures for no ejection. 
The results reported here show that, associated with these reduced pressures, 
there are greatly reduced values of heat flux. However, in the vicinity of the 
stagnation circle the pressures are somewhat above their values for no ejection. 
Associated with this there is an increased heat flux in a region around the stagna- 
tion circle. The nett result is that, although the region of severest heating is 
shifted, the overall heating is not greatly reduced, at least not unless large 
coolant flow rates are employed. 

The second flow régime occurs when the coolant flow rate is sufficiently small 
that the coolant gas does not cause the bow shock wave to bulge out. With the 
external flow sensibly undisturbed, relatively large reductions in heat flux are 
obtained. The reductions in total heat transfer over the body are of the order 
of the ‘heat capacity’ of the coolant gas, defined for this purpose as the amount 
of heat that the gas can absorb in having its temperature raised from that just 
before ejection to the stagnation temperature of the main air stream. 

The work reported here is an experimental determination of the local rates 
of heat transfer over a bluff body for various flow rates of a coolant gas, either 
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nitrogen or helium, ejected either straight out or with swirl, in the manner 
already described. The momentum flow rate of the ejected coolant was not 
varied independently of the mass flow rate. To doso would be one way of extending 
the present work, for it would be extremely interesting to eject coolant gas at 
larger rates of mass flow, but with insufficient momentum to cause the bow 
shock wave to bulge out. 


2. Brief description of the experiments 


The experiments were made in the 5 in. x 5 in. hypersonic wind tunnel of the 
Guggenheim Aeronautical Laboratory at the California Institute of Technology. 
The nominal Mach number was 5-8, and the Reynolds number per inch based 
on free-stream conditions was about 0-2 million. At this Reynolds number and 
Mach number it was assured that the boundary layer on the model was laminar. 
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FicuRE 1. Construction of models. (a) Pressure model. (b) Temperature model. 


The shape chosen for investigation consisted of a portion of a cone of 10° 
semi-angle, with the apex rounded to form a portion of an oblate spheroid of 
eccentricity 0-7315. Two models were constructed: one to obtain pressure 
distributions and the other to obtain temperature distributions. Figure 1 shows 
the internal construction of the models. 
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The pressure model was constructed of brass with nine conventional static 
pressure orifices of 0-016 in. diameter. Figure 2 shows the position of the 
pressure orifices, which were located helically on eight meridians. All pressures 
have been quoted as pressure coefficients, Cp, defined by 


where P,, p,, and V,, are respectively the pressure, density and speed of the 
freestream. Cp has been normalized by dividing it by the pressure coefficient, 
Cp,, at the nose for zero incidence and with no coolant ejection ; which is, of course, 
the pressure coefficient behind a normal shock at the tunnel Mach number. 






__| All dimensions 
“1J_}_ are in inches 








Relative location of Relative location of 
pressure orifices thermocouples 


FIGURE 2. Geometry of models. 


The temperature model was constructed of a mild steel with a wall thickness 
of 0-050 in., and coated with a layer of porcelain of thickness 0-020 in. Along 
one meridian of the temperature model nine thermocouples were formed by 
cementing a constantan wire of 0-001 in. diameter on the model, and then 
removing the cement coating at nine precisely located points with a knife 
blade. It thus became possible to make contact with the wire by means of silver 
paint at these nine points, thereby forming nine silver-constantan junctions 
(see figure 2). Painted silver leads connected the thermocouples so formed to 
turret posts set in fibre-glass on the base of the model. Some heat-meters also 
were installed on the temperature model with which it was hoped to measure 
the heat flux directly. However, unsatisfactory results were obtained with 
them, which will not be reported here. 

The temperature mode! was cooled by water, the inflow being within the 
return-flow (see figure 1). The flow rate was on the average about 0-17 U.S. 
gallons per minute, varying at most by a factor of 2. The inflow water tem- 
perature whenever measured was always within | °F of 76 °F, and the maximum 
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measured temperature rise overall, that is due to the heat picked up not only 
in the model but also in the attendant piping, was never greater than 7 °F, 
and usually less than 3 °F. It can be shown (Warren 1958) that the heat flux, 
g, at a particular point across a thin skin of low thermal conductivity is fairly 
uniform, and is given very closely by 

= (T— 7). 


q = 


where / is the thermal conductivity of the skin material, 6 is the thickness of 
the skin, 7’ is the temperature of the outer surface at the particular point, and 
T, is the corresponding temperature of the inner surface. 

Since the thermal conductivity of porcelain is about one-sixtieth of that of 
steel, it is not unreasonable to assume that most of the temperature drop 
across the wall of the model occurs in the 0-020 in. layer of porcelain, and 
that the inner surface of the porcelain is not far removed from the water tem- 
perature, which, in the model, was about 80 °F. In view of the failure to obtain 
satisfactory results with the heat-meters, normalized local rates of heat transfer 
were obtained from the measurements of surface temperature by assuming that 


qj | —_ 7. 


i 

do T%-T; 
where g, is the heat flux at the nose for zero incidence and with no coolant 
ejection, 7’ is the measured surface temperature, 7, is the surface temperature 
at the nose for zero incidence and with no coolant ejection (measured to be 
130 °F), and 7, is the inner surface temperature, assumed to be 80 °F. 

The nitrogen and helium gases used in the ejection studies were fed into the 
model via a pipe within the inflow water pipe. This procedure ensured that the 
coolant gas reached the model at a fairly constant temperature. Within the 
model the gas entered a plenum chamber (see figure 1) where a pressure probe 
and an iron-constantan thermocouple enabled the pressure and temperature 
of the gas to be determined. From the plenum chamber the gas was conveyed 
to the nose of the model through an ejection pipe of 0-081 in. inside diameter 
and 0-5 in. long, the edges of the ejection hole being rounded to a radius of 
about 0-02 in. The same physical parts for the coolant ejection were used in both 
the pressure and temperature models. The coolant mass flow rate, 7, has been 
quoted as a mass flow coefficient, C,,, based on the free-stream mass flow rate 
through a capture area equivalent to the frontal area of the model. so that 


where D is the model base diameter. 
The swirl that was given to the coolant in some of the tests was obtained by 
inserting in the ejection pipe a twisted strip of brass 0-01 in. thick and 0-5 in. 


long. The strip was twisted through one complete turn clockwise, so that the 
helical advance angle at the wall was about 63°. It was appreciated that such 
a device would not produce initially a rotationally symmetric coolant flow, 
but it was hoped that, once away from the ejection orifice, the flow would 
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even out. On the pressure model the setting of the swirler at the ejection orifice 
made an angle of about 40° with the meridian through pressure orifices | and 9, 
as indicated in figure 2. On the temperature model two positions were tested: 
one made an angle of about 45° with the meridian through the thermocouples, 
as indicated in figure 2, and the other made an angle of about 90°. 

Further details of the experiments, together with all the results obtained, are 
given elsewhere (Warren 1958). Only a selection of the results are presented 


here. 


3. Results for no ejection 

Distributions of pressure and heat flux in the plane of incidence with no 
ejection at incidences of 0, 4 and 8° are shown in figures 3 and 4, in terms of the 
distance, s, from the model nose along a meridian. It was found that there is 
no noticeable variation of pressure or heat flux with incidence in the meridian 
plane normal to the plane of incidence, at least up to an incidence of 8°. 
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FicurE 3. Distributions of pressure in the plane of incidence for no ejection (C;, = 0). 
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The pressure results are in fair agreement with Newtonian theory, which 
states that the pressure is given by 
Y 

Cp 


( = sin* 7), 


Po 
where 7 is the angle of the surface to the direction of motion, except on the 
conical portion, where the Newtonian approach underestimates the pressure. 
The pressure results are also in good agreement with those of Richards (1957), 
who earlier tested a model of nominally identical geometry. 

The heat flux results at zero incidence are in fair agreement with calculations 
based on the theory of Lees (1956). In these calculations the measured, as 
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distinct from the theoretical, pressure distribution was used, and allowance 
was made for the small effects associated with a pressure gradient. The poor 
agreement over the conical portion may be associated with the measurements, 
for the method whereby the heat flux is derived is of least accuracy when the 
heat flux itself is small. Calculations based on the theory of Cohen & Reshotko 
(1956), although more difficult to perform, were in substantial agreement with 
those based on the theory of Lees. 
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Figure 4. Distributions of heat flux in the plane of incidence for no ejection (C,, = 0). 
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4. Variation of ejection conditions with coolant flow rate 


On the assumption of one-dimensional flow in the coolant gas pipe, the 
Mach number in the plenum chamber never exceeds 0-03, based on the measure- 
ments of mass flow rate and of pressure and temperature in the plenum chamber. 
The plenum chamber pressures and temperatures can therefore be treated as 
representing total conditions. On this assumption, and again that of one- 
dimensional flow, the Mach number, /,, at the entry to the final ejection pipe 
can be calculated. This Mach number is shown in figures 5 and 6. The constancy 
of Mach number shown in figure 5 for mass flow coefficients above about 0-005 
for nitrogen or 0-0015 for helium is evidence that, for ejection with swirl, the 
flow in the swirl pipe is choked. Figure 6 shows that for straight-out ejection 
the onset of choking is delayed to a mass flow coefficient about 0-008 for 
nitrogen and 0-0025 for helium. The Mach numbers themselves should be treated 
as qualitative only, owing to the assumption of one-dimensional flow made in 
their derivation, and to the fact that Mach number is extremely sensitive to 
area ratio near a Mach number of unity. 

For straight-out ejection the momentum flow rate of the ejected gas can be 
calculated, again on a one-dimensional flow assumption, and ignoring any 


losses in the final ejection pipe, which are estimated to be small. The calculated 
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momentum flow rates are shown in figure 7 as momentum flow coefficients 
C.,, defined as ee 
B 0. = mV,; 

"Po Vim(2D)’ 

where J, is the velocity of ejection. It will be noted that for small mass flow 
rates, such that the flow in the ejection pipe is not choked, the momentum 
flow rate for helium, and hence its speed of ejection, is about seven times that 
for nitrogen for the same mass flow rate, as would be expected from the relative 
densities of helium and nitrogen. 
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FIGURE 5 FIGURE 6 
FicurE 5. Mach number of coolant at entry to the final swirl pipe for ejection with swirl. 
0 @.a=0°:A 4,a=4°; Om, «a=8°. The flagged points are for the pressure model, the 


unflagged points are for the temperature model. 


Ficure 6. Mach number of coolant at entry to the final ejection pipe for straight-out 
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FicurE 7. Variation of coolant momentum flow with coolant mass flow for straight-out 


ejection. O @,a=0°;A 4,42=4;08,a=8. 





408 C. H. E. Warren 


5. Effects of ejection on the flow pattern 


The effects of ejection on the flow pattern can be deduced from an examina- 
tion of the schlieren photographs and the measured surface pressure dis- 
tributions. 

For straight-out ejection at very small flow rates no noticeable change in 
the schlieren picture occurs (see figures 8a,b, plate 1). At a certain critical flow 
rate, however, the bow shock wave begins to bulge out (see figure 8c, plate 1). 
This critical flow rate appears to increase with increase in incidence (see 


figures 84,h, plate 2). At zero incidence it occurs at a mass flow coefficient of 


about 0-003 with nitrogen and of about 0-001 with helium. From figure 7 we 
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FicurE 9. Distributions of pressure at zero incidence (« = 0°) for straight-out ejection. 
(MeMahon’s results for a sphere-cone model.) C,, =m/p,V.7R?. 


see that these mass flow coefficients correspond to a momentum flow coefficient 
in each case of about 0-00035,* thereby confirming the work of Lam (1959) 
who suggests that the flow pattern is determined, at small flow rates, by the 
momentum flow coefficient of the ejected gas. Above the critical flow rate the 
bow shock wave is bulged out, and from the results of McMahon (1958) there 
is an associated marked change in the surface pressure distribution. McMahon 
studied the same problem on a model of about the same size, but having a 
spherical instead of a spheroidal nose, and an ejection orifice smaller by about 
25°,. His tests were made at the same Mach number and about the same 
Reynolds number. His results are shown in figure 9. It will be noted that the 
surface-pressure distribution is not markedly changed for a mass flow coefficient 
of 0-002 with nitrogen, corresponding to an unbulged bow shock wave. For 
all other mass flow coefficients shown, however, there are marked changes, 
corresponding to bulged bow shock waves. The pressure distributions are in 
fact similar to those that occur when a spike protrudes from the nose of a model, 


* This momentum flow coefficient based on model base area is equivalent to a momen- 


tum flow coefficient of 0-12 based on ejection orifice area. 
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implying that the ejected gas acts as a ‘gas spike’. It will be seen from figure 9 
that ejection reduces the pressures in the vicinity of the nose. At the smaller 
flow rates there is, somewhat surprisingly, an associated reduction in the 
plenum chamber pressure. P,, (see figure 10). Moreover, the plenum chamber 
pressure varies appreciably with incidence at the smaller flow rates. This is 
presumably because the conditions at the ejection orifice vary appreciably with 
incidence. With helium when the mass flow coefficient is above about 0-003 
this dependence on incidence disappears: this is no doubt associated with the 
fact that above this flow rate the flow in the ejection pipe is choked, as men- 
tioned in §4. With nitrogen the plenum chamber pressure depends on incidence 
at all flow rates up to that corresponding to a mass flow coefficient of 0-008, the 
maximum tested: this too is in keeping with the result, obtained in §4, that the 
flow in the ejection pipe is unchoked for mass flow coefficients less than 0-008. 
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FicurE 10. Plenum chamber pressure for straight-out ejection. 
0@,a=0; AAadA,a=*;s O8,a4=S8. 


For ejection with swirl the critical flow rate above which the bow shock 
wave begins to bulge out is greater than for straight-out ejection (see figures 
8c,d,i and j, plates 1, 2). It would appear, therefore, that the swirler is achieving 
its desired effect. At zero incidence the critical flow rate occurs at a mass flow 
coefficient of about 0-008 with nitrogen and of about 0-003 with helium. For 
ejection with swirl the concept of a momentum flow rate is not so easy to 
formulate, but it is interesting to note that the ratio of the two critical mass 
flow rates for ejection with swirl is roughly the same as the corresponding ratio 
for straight-out ejection, namely about 7:1. For ejection with swirl, however, 
it was found that the critical flow rate decreases with increase in incidence 
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(see figures 8d, e and f, plate 1). The variations in the surface pressure distribution 
for ejection with swirl (as measured in the model of the present investigation) 
are shown in figures 11 and 12. As noted also in the results for straight-out 
ejection the surface-pressure distributions at zero incidence (figure 11) are not 
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FicurE 11. Distributions of pressure at zero incidence (« = 0°) for ejection with swirl. 
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FicurE 12. Distributions of pressure in the plane of incidence 
for the ejection with swirl of nitrogen. 


markedly changed for mass flow coefficients for which the schlieren pictures 
indicate an unbulged bow shock wave (mass flow coefficients less than about 
0-008 with nitrogen and less than about 0-003 with helium), but there are 
marked changes corresponding to bulged bow shock waves. 

The variations in plenum chamber pressure are shown in figure 13. The 
plenum chamber pressures do not vary noticeably with incidence: this is no 


doubt because, with swirl, the flow in the swirl pipe is choked when the mass 
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flow coefficient is above about 0-005 for nitrogen or 0-0015 for helium (see § 4); 
and up to these flow rates, as we have seen, the bow shock wave is not bulged 
out, the flow ahead of the body is little affected, and the pressures in the 
vicinity of the nose do not vary much with either incidence or coolant flow rate. 
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FIGURE 13. Plenum chamber pressure for ejection with swirl. © @,«=0°; A 4,a=4°; 
|] g@, a= 8°. The flagged points are for the pressure model, which was uncooled; the 
unflagged points are for the temperature model. 


It will be noted, however, by comparing figures 10 and 13 that, for the same 
mass flow rate, much greater driving pressures are required for ejection with 
swirl than for straight-out ejection. The slightly greater driving pressures 
required with the pressure model compared with the temperature model, 
indicated in figure 13, are because the pressure model was uncooled, and 
accordingly the temperatures in the plenum chamber were some 100 °F higher 


than with the temperature model. 
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6. Effect of ejection on the distribution of heat flux 

The distributions of heat flux for ejection with swirl at zero incidence are 
shown in figure 14. Results are shown for one flow rate with each coolant, the 
flow rates chosen being ones for which the schlieren and pressure studies con- 
sidered in §5 indicated that the swirler was working, in the sense that the bow 
shock wave was not bulged out (see figures 8d, j, plates 1, 2) and the distributions 
of pressure were not markedly different from the case of no ejection (see 
figure 11). The first thing that one notices in figure 14 is that the results depend 
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FicurE 14. Distributions of heat flux at zero incidence (a = 0°) for ejection with swirl. 
Nitrogen: ©,C,;, = 0; A, C,, = 0-006, swirler at 45°; 7, C,, = 0-006, swirler at 90°. 
Helium: ©, C;,=0; A, C,, = 0:002, swirler at 45°; V7, C,, = 0-002, swirler at 90°. 


critically upon the setting of the swirler, implying that the swirler is not giving 
a rotationally symmetric distribution of coolant gas over the model. In view 
of these results an attempt should have been made to obtain a more nearly 
rotationally symmetric distribution of coolant gas. However, the main obser- 
vation obtained from figure 14, and from other results not shown here, is that 
ejection with swirl seems to lead to an increase in heat flux. Accordingly the 
matter was not pursued, but attention directed instead to straight-out ejection. 

An explanation of the increase in heat flux produced by ejection with swirl 
is provided by the Reynolds analogy, which states that the ratio of the Nusselt 
number to the product of the surface friction coefficient and the Reynolds num- 
ber is a constant depending upon the surface temperature and the longitudinal 


yressure gradient, or Pin 2 pies 
] 5 k(7 _ T') 


Lh, 


where 7’ is the surface temperature, 7), is the adiabatic surface temperature 
(that is the surface temperature for which the heat flux is zero), I, is the 
velocity at the external edge of the boundary layer, 7 is the surface shear stress, 
k is the thermal conductivity, and uw is the dynamic viscosity. 
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Now the main effect of the introduction of a coolant is to reduce the driving 
temperature difference (7',—7'), and this is the basis on which it works, but in 
ejection with swirl the coolant is injected tangentially into the boundary layer 
at the nose with very high velocity, much greater than the average velocity of 
the uninjected boundary layer, with the result that the surface shear stress is 
increased manyfold. The results obtained seem to indicate that the increase 
in 7 outweighs the reduction in (7',— 7') and the changes in the other quantities, 
particularly * and y. If tangential injection into a boundary layer is to be 
instrumental in reducing heat flux it must be done with much less momentum 
for a given mass flow than was achieved in these experiments. 
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FicureE 15. Distributions of heat flux at zero incidence (a = 0°) 
for straight-out ejection. 


The distribution of heat flux for straight-out ejection at zero incidence are 
shown in figure 15. The curve for nitrogen at a mass flow coefficient of 0-002 
is of a markedly different shape from the curves for the other three cases. 
Moreover, it corresponds to a reduction in heat flux along the entire body profile. 
From the schlieren and pressure studies considered in §5 this case is one for 
which the coolant is introduced with sufficiently little momentum that it does 
not cause the bow shock wave to bulge out (see figure 8b, plate 1), or the pres- 
sure distribution to be markedly different from the case of no ejection (see 
figure 9). The other three cases (nitrogen at a mass flow coefficient of 0-006 and 
helium at mass flow coefficients of 0-002 and 0-004) all correspond to bulged 
bow shock waves (see figures 8c,i, plates 1, 2), and markedly different pressure 
distributions (see figure 9). In these cases although large reductions in heat 
flux are obtained in the vicinity of the nose, there is an increase further round 
the profile, where s/D is about 0-25. These results agree with those of McMahon 
(1958). He found large reductions in the average heat flux over the nose sector 
of a sphere of 60° included angle, but considerably smaller reductions over a 
sector of 120° included angle. 
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The explanation of these results is provided by considering the associated 
flow patterns. As shown in §5, in these cases the ejected gas acts as a ‘gas spike’, 
causing the main air flow to separate ahead of the model thereby leaving a 
roughly conical ‘dead-air’ region in the vicinity of the nose. The flow external to 
the boundary layer, instead of having a stagnation point at the nose, now has 
a ‘stagnation circle’ where the incident flow meets the body, along which the 
heat flux is increased compared with the heat flux at the same station for no 
ejection. Inside the dead-air region the pressure is greatly reduced (see figure 9), 
and this reduction in pressure is largely responsible for the reduced heat flux 
in the vicinity of the nose. 

It is of interest to consider the ratio of the reductions in total heat transfer 
to the ‘heat capacity’ of the ejected gas, taken as 

met (i,—T7,)- 


p 


This expression represents the amount of heat that the gas can absorb when its 
temperature is raised from that in the plenum chamber, 7,,, to the stagnation 
temperature of the main air flow, 7,; c, is the specific heat of the gas at 
constant pressure. For nitrogen ejected straight out at a mass flow coefficient 
of 0-002, corresponding to an unbulged bow shock wave, the ratio of the 
reduction in total heat transfer up to the station s/D=1 to the ‘heat 
capacity’ of the ejected gas is about 1-1. That this is greater than unity 
implies that the nitrogen is acting more than as a blanket; that is, a greater 
reduction in heat transfer is being achieved than would correspond to the 
absorption by the nitrogen of the heat that would otherwise have been trans- 
ferred to the model. The ratio of the reduction in total heat transfer up to the 
shoulder (s/D = 0-459) to the ‘heat capacity’ of the ejected gas is about 0-4, 
which is still remarkably high. These figures should be compared with those for 
helium for the same mass flow rate, remembering that in this case the shock 
wave is now bulged out (see figure 87, plate 2). The ratio of the reduction in 
total heat transfer up to the station s/D=1 to the ‘heat capacity’ of the 
ejected gas is about 0-1, which is less than one-tenth of the value with nitrogen: 
the ratio based on the total heat transfer up to the shoulder is roughly zero. 
The implication of these results seems to be that if the gas can be ejected 
without causing the bow shock wave to bulge out, as is the case with nitrogen 
ejected at a mass flow coefficient of 0-002, then a steady reduction in heat 
flux with increase in flow rate should be achieved. This state of affairs can 
be achieved by arranging that the mass flow rate be increased without in- 
creasing the momentum flow rate, so that on impingement with the air the 
coolant is directed around the body without too much mixing and consequent 
reduction of its heat alleviation properties. The importance of this observation 
is brought out by the results for nitrogen and helium ejected at a mass flow 
coefficient of 0-002 already referred to. Although it has five times the heat 
sapacity and is one-seventh of the density, implying a thicker layer for the 
same mass flow rate and velocity, helium is not as effective as nitrogen in re- 
ducing the heat flux at this flow rate simply because its heat alleviation 
properties are enormously reduced by mixing. It would obviously be extremely 
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interesting to extend these experiments by varying the ejection pipe in such 
a way that the mass flow rate with a given gas could be varied independently 
of the momentum flow rate. In this way it should be possible to delay bulging 
the bow shock wave to higher mass flow rates, and substantial reductions in 
heat flux for quite small mass flow rates should thereby be achieved. 

The effect of incidence on the distribution of heat flux for straight-out 
ejection is shown in figure 16 for nitrogen. Qualitatively the results are similar 
to those at zero incidence. At a mass flow coefficient of 0-002, corresponding to 
an unbulged bow shock wave, the heat flux is reduced less on the windward 
meridian than on the leeward meridian; and at an incidence of 8° ejection is 
having very little effect on the windward meridian. At a mass flow coefficient 
of 0-006, corresponding to a bulged bow shock wave, the heat flux is greater, 
over most of the model, than the heat flux at the lower flow rate, and on the 
windward meridian it is greater than with no ejection. 
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Figure 16. Distributions of heat flux in the plane of incidence 
for the straight-out ejection of nitrogen. 


7. Conclusions 

In the course of the investigation three ways of ejecting a coolant gas from 
the nose of a bluff body have become apparent, and something has been learned 
about the heat transfer rates associated with them. 

The method that was tried initially was to encourage the coolant gas to flow 
tangentially to the model surface at ejection by giving it some swirl in the 
ejection pipe. The method was successful in the sense that, at the smaller flow 
rates, the coolant gas was ejected with very little disturbance to the flow 
external to the boundary layer, as witnessed by the pressure distributions and 
schlieren studies. However, this method of ejection is not satisfactory as a 
means of easing the heating problem. Although the cool layers of coolant gas 
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reduce the driving temperature difference the high tangential velocity with 
which the coolant is injected into the boundary layer increases the surface shear 
considerably, and hence, by the Reynolds analogy, the heat flux. The results 
indicate that the effect of the increased shear predominates, at least for the 
conditions considered here. This method of ejection was therefore abandoned. 

Attention was then centred on straight-out ejection, but here two flow 
régimes were found. Except at the smaller flow rates the coolant gas acts 
as a ‘gas spike’, causing the bow shock wave to bulge out and the main air flow 
to separate ahead of the model, thereby leaving a roughly conical ‘dead-air’ 
region in the vicinity of the nose. The flow external to the boundary layer, 
instead of having a stagnation point at the nose, now has a ‘stagnation circle’. 
Inside the dead-air region the pressure is greatly reduced compared with the 
case of no ejection, and, associated with this reduced pressure, there is a greatly 
reduced heat flux. However, in the vicinity of the stagnation circle the pressure 
is somewhat above its value for the case of no ejection, and associated with this 
there is an increased heat flux. The nett result is that, although the region of 
severest heating is shifted, the overall heat transfer rate is not greatly reduced, 
at least not unless large flow rates are employed. 

At the smaller flow rates the momentum of the ejected gas is not sufficient 
to disturb the external flow appreciably. The critical flow rate corresponds to 
a momentum flow coefficient of about 0-00035. With nitrogen this corresponds 
to a mass flow coefficient of 0-003, and with helium to a mass flow coefficient of 
0-001. This method of ejection was investigated with nitrogen at a mass flow 
coefficient of 0-002, and under these conditions relatively large reductions in 
heat flux were obtained all over the model. The reductions in total heat transfer 
were of the order of the ‘heat capacity’ of the coolant nitrogen, where the ‘heat 
capacity’ is defined for this purpose as the amount of heat that the gas can 
absorb in having its temperature raised from that just before ejection (i.e. in 
the plenum chamber) to the stagnation temperature of the main air stream. In 


comparison, for the straight-out ejection of helium at the same mass flow 


coefficient of 0-002, but a case in which the bow shock wave is now bulged out, 
the reduction in total heat transfer over the model was no more than one-tenth 
of the ‘heat capacity’ of the coolant helium. 

Clearly the method of ejecting straight out, but without sufficient momentum 
to cause the blow shock wave to bulge out, seems extremely promising as a 
means of alleviating the heating problem. It would form an obvious and 
interesting extension of the present experiments to be able to vary the momen- 
tum flow rates and mass flow rates of an ejected coolant independently. This 
could be done by varying the size of the ejection pipe. 


This investigation was made under Army Ordnance Contract No. DA-04-495- 
Ord-19 whilst the author was at the Guggenheim Aeronautical Laboratory, 
California Institute of Technology. on a Harkness Fellowship of the Common- 
wealth Fund of New York. 
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(j) 
Figure S (plate 2). Sehlieren photographs. (g) Nitrogen ejected straight out, C. = 0-004; 
| | | { J ee 


x=. (h) Nitrogen ejected straight out, C, = 0-006, 2=8. (4) Helium ejected straight 
out, C. = 0-002, a2 = 0°. (yj) Helium ejected with swirl, C.. = 0-002; a = 0°. 
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On Kelvin’s ship-wave pattern 


By F. URSELL 


Department of Applied Mathematics and Theoretical Physics, University of Cambridge 
(Received 25 October 1959) 


When a concentrated pressure travels with constant velocity over the free surface 
of water, it carries with it a familiar pattern of ship waves. Let viscosity and 
surface tension be neglected, let the free-surface condition be linearized, let the 
depth of water be assumed infinite, and let initial transient effects be ignored. 
Then, as is well known, the wave motion everywhere can be found by standard 
methods in the form of a double integral. The wave pattern at a great distance 
behind the disturbance can be found by an application of the ordinary method 
of stationary phase, which shows that the wave amplitude is considerable 
inside an angle bounded by the two horizontalrays@ = +40, from the disturbance, 
where (7, = sin-!} = 194°. But the method fails in two regions, near the track 
&# = 0 of the pressure point, and near the critical lines 0 = + 8,. 

These two regions are treated in the present paper. It is shown that near 
? = 0 the linearized surface elevation oscillates with indefinitely increasing 
amplitude and indefinitely decreasing wavelength. (This result holds only when 
the pressure is concentrated at a point and applied at the free surface.) Near 
the critical lines the surface elevation at a greater distance behind the pressure 
point can be expressed in terms of Airy functions, and this expression goes over 
into the known wave pattern inside the critical angle. It is shown that near the 
critical lines the crest length increases as the cube root of the distance, and that 
the separation between crests remains constant. Contour maps of the wave sur- 
face are given for three distances behind the moving pressure point. 


1. Introduction 

When a disturbance (e.g. a ship) travels on a water surface, it carries with it a 
familiar pattern of bow and stern waves which was first explained mathematically 
by Lord Kelvin (Sir W. Thomson 1891). Instead of ship waves he considered 
the waves generated by a prescribed pressure distribution moving with a constant 
velocity U and acting on the water surface. The magnitude of the pressure and 
the resulting wave slope were assumed to be so small that the equations of motion 
could be linearized; viscosity and surface tension were neglected. Under these 
assumptions it is sufficient (in principle) to calculate the waves due to a moving 
concentrated pressure point; the effect of a distribution of surface pressure can 
thence be found by integration. Kelvin showed, apparently by using his principle 
of stationary phase, that the characteristic wave pattern of a pressure point is 
the superposition of two sets of waves effectively confined between the two 


vertical planes ) = +6, = +194, where the co-ordinates are taken as in § 2 below; 
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for a figure see Lamb (1932, §256). A simple argument using such ideas was 
given by Havelock (1908, p. 417), and by Lamb (1932, §256). All these results 
are valid at a distance of many wavelengths behind the pressure point; near the 
pressure point the principle of stationary phase is not applicable. It was found 
that near the critical lines 0 = + 0, the amplitude is much larger than elsewhere 
at the same distance from the pressure point, and the simplest form of the prin- 
ciple of the stationary phase could not be used there. (Kelvin predicted infinite 
amplitudes on 7 = @,; the amplitude variation on 0 = @, was first given correctly 
by Havelock (1908, see equation (4.1) below).) The study of the wave pattern 
near, but not actually on, a critical line is more difficult and was undertaken by 
Hogner (1923) who was able to write the surface elevation as the sum of a double 
and of a single integral, and to show that only the latter is significant at a large 
distance behind the pressure point. This is convenient because the application 
of the methods of steepest descents and of stationary phase to single integrals is 
simple. There are, however, some difficulties in his work. First, the single integral 
has one finite limit of integration, and this leads to complications (see p. 19 of 
his paper). Secondly, near the lines 0 = + @,, where a single asymptotic expres- 
sion in terms of Airy functions would be expected (cf. equation (4.12) below), 
Hogner’s expression changes its form as @ passes through 0, (Hogner, 1923, 
equations (82), (83)). 

A new treatment of the problem has been given by Peters (1949). The surface 
elevation was obtained by him in a form similar to Hogner’s, as the sum of a 
double and a single integral where again only the latter contributes significantly 
to the waves far behind the pressure point; but Peters’s integral, unlike Hogner’s, 
has both limits of integration at infinity. By an immediate application of the 
method of steepest descents he found expressions which are valid when @ is 
not near 9,. An expression was also given for 7 on and near 6 = @,; but this is 
valid only when N?|@—0,| is small, see §4 below. 

While the work of Peters has thus reduced the problem of finding the waves 
far behind the pressure point to the asymptotic evaluation of a single integral 
containing a large parameter, there are two regions which have not yet received 
adequate treatment, and these will be considered in the present paper. One is 
the vicinity of the track of the pressure point where @ and z are small. The limits 
are non-uniform: we shall see by applying the method of steepest descents that 
the amplitude tends to oo as @ tends to 0 on z = 0 (but on 0? = 0 it is finite as 
z tends to 0). The other region is the neighbourhood of 6 = 0,, where an expression 
in terms of Airy functions will be given; here the method of Chester, Friedman & 
Ursell (1957) is applicable, and the expression is valid in some finite angle 
including the line 0 = @,. 


2. The elevation due to a travelling pressure point 

Peters has calculated the velocity potential and the surface elevation 
C(x, y; U,g) above the mean water level. In this calculation the exact free- 
surface condition has been replaced by its linearized approximation. Take the 
origin of co-ordinates travelling with the surface pressure point, the x-axis 
horizontal along the track of the disturbance, and the y-axis horizontal at 
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right angles to the z-axis. Write AK =g/U? and define N and @ by the 


equations ; . . — 
: Kx =Ncos0, Ky = Nsin0. (2.1) 
The critical angle @, is defined by 
as ok, ae 
6, = sin} = tan-!5—5 = 194°; 
2/2 


and the total force acting on the fluid is denoted by Py. 
Then Peters’s expressions for €(x, y) are: 
for x > 0 


C(x, y) = G(2, y) + &.(2, y), 
where 
wal? 00 , wf sin’ y cos (ky cos y) 
C(x, y) = dk k? e-kx d) i =.» 
Po seit J0 J0 K? + k? sin? y 

and 

1PgY 
a => Salt, 4 

P, K2o2 »Y) 


y* 
z->0¢ 


= 2lim | (1+?) cos {Ayu (1+ u?)} sin {Ka , (1+ u?)}exp {— K2(1+u?)} du 
0 
Px 
=im lim | (1+ u?)exp [iV {(cos 0 — usin) .(1+u?)!}exp {— K2(1+ wu?) du; 


200 L (2.2) 


while for x < 0, (x,y) = €(-—2,y). At the pressure point x = 0, y = 0 the inte- 
grals diverge. 

It appears impossible to express ¢, and ¢, in terms of known functions, but 
asymptotic expansions for the amplitude at a large distance (N > 1) behind 
the pressure point can be found, and these show that the contribution of ¢, is 
negligible inside the critical angle || < 0,. For according to P, p. 142,* 


9 roy’o . 9 . al 
7? pK?U2 | 7 37y” . _ [27 , cos (ky cos 
"P Cetra © es _ § “4 ” He | dy = 2. 
0 (a°+y")? 2(r*+y")? Jo 0 a a 
(2.3) 
and when the substitution k = Ausiny is made in the double integral, this be- 
comes 
Por ehn (K j 3 7,4 
2 K sin y)? u es ee 
| du | dy oe . ) exp (— Axusin y) cos (Kyusiny cos y) 
-u? 
/0 70 : 
a iy Oe ut ee ee 
= = | dy(Ksiny)? | du baye exP (—Kausiny —iKyusiny cos y) 
2 +u* 
J0 Jo 
Lor? 
= >| dy(K siny)? F)(Kasiny +ikysin y cos y), (2.4) 
= 0 
where by definition 
*D 9,4 p—uZ » C2 p—uZ 
ure 2. J e 
F,(Z) = du = —.—=+ du. 
. Jo 1+u? PE je 1+ 


* References preceded by the letters H or P are to Hogner (1923) or Peters (1949), 
respectively. 
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The last integral is related to the exponential integral (Jeffreys & Jeffreys 1946, 
p. 443), and it can be shown that for re Z > 0 and |Z| < 1 we have 
| Fo(Z)| < A|Z|-3; 
while for re Z > 0 and |Z| > 1 we have |F,(Z)| < A|Z|->, by the method of 
steepest descents. To find a bound for (2.4), write it as 





a5 


ig if ifs 
| + +5] = I,+1,+I , say. 


9 9 of. 
67 


-/0 ~~/6 


In J, introduce the new variable 6 = Kxsiny+iKysiny cosy; then 














iy = |Kxcosy+ikycos2y| > 4|Kx+iky| = 3N, 
and |Ksiny| = |A(x+iycosy)“| < 28 |a+iy|-! = 22KN-. 
TT ‘ ’ — 
Thus, [Zi| <3 |28|8K3N-3| F,(f)|2N—-|dp| < AK3N-4. (2.5) 
ae 
Similarly, |I,| < AK®N-*. 


To find a bound for J,, deform the path of integration in the y-plane to pass 
below y = 47; on the new path, the real part of Kxsiny+iKysiny cosy is 


positive, and > A|Kr+iky| 








|Kxsin y +iKysin y cos y 


which is uniformly large when JN is large. Thus 


lT, < AK®N-5, (2.6) 








<A [ ~ K8N-5|dy 
and equations (2.5) and (2.6) show that the double integral in (2.3) is O(K3N~). 
Thus ¢, is of order N-* everywhere, while ¢, is of order V-? inside the critical angle, 
and exponentially small outside the critical angle, see P, equations (5.2), (5.4). 
Thus, ¢, is dominant inside the critical angle, ¢, outside. We are considering the 
asymptotic form of ¢(2,y) for large N = K|x+iy|; the expressions (2.3) and 
P (5.2), (5.4) are adequate except near 6 = 0 and 6=+40,. (They will not be 
quoted here; see, however, §4 below.) We shall now obtain new expansions valid 
near 6 = 0 and @ = 6,. 


3. The track of the pressure point 

We shall have to evaluate (2.2) for small 7, where, by reason of the symmetry 
of the pattern, we shall take 6 > 0. If we attempt to find the limit in (2.2) by 
putting z = 0 in (2.2) we see that the resulting integral does not converge. This 
difficulty is avoided if the path of integration is deformed before z is placed equal 


to 0. The new path is chosen so that 

exp [tN {(cos @— usin @) /(1+u?)}] (3.1) 
tends to zero rapidly as |u| — oo along the new path, and there are two cases, 
# = 0 and @ > 0, which must be considered separately. 


(i) If 6=0, take a path of integration along the straight line from 
u=—ooexp(imi) through w=0 to cexp(4zi), say. Along this path 


> 
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exp {iN /(1+w?)} > 0 as |u| > », and exp{—K2(1+4?)} is uniformly bounded 
for all w on the path and all z > 0. Thus, this integral is uniformly convergent as 
z-> 0, and 


ex exp(}z7i) 
(x,0) = im (1+u*)exp {iN /(1+u*)} du. (3.2) 


1 
. Lexplgrd) 


T™PJ . 
8 ry Se 
P, K* 
(Clearly a good deal of latitude is allowed in the choice of the path.) This may be 
written in terms of Bessel functions, (P, p. 140, ef. Hogner 1924, p. 7): 
PM oe ce a P,K2( 2 \3 
. £0" 3Y (Kx) —Y,(Kx)} ~ —-° 
4pg , pg \wke 
and so €, is finite when Kz is positive. 
(ii) If0>0, take a path of integration along the straight line from 
u = —wexp(imi) tou = 0, and thence along the straight line to uw = cexp(— §772) 
Along this path the expression (3.1) tends to 0 as |w| +00, and exp {— Kz(1+ u*)} 
is uniformly bounded. Thus, as in (i), the integral is uniformly convergent as 


cos(Ka— 47), (3.3) 


z > 0 for fixed 0, and 


7pq . : P20 exp(—sz7t) P 
~ ~5 Go(2, y) = 1M (1 +wu*)exp 
Pon* 





iN{(cos @—wsin #) /(1+u*)}] du. 
(3.4) 
We note that the limits of integration in (3.2) and (3.4) are different, and we are 
led to suspect singular behaviour near # = 0. This is confirmed by the following 
calculation. First, deform the path of (3.4) into two paths of steepest descent 
through the two saddle points. These points are the roots of 
(d/du) {(cos 9 — usin A) (1 +u?)} = 0, 
i.e. the roots of u cos @—(1+2u?)sind = 0, 
i.e. the points 
u,(A) = }{eot A+ ,/(cot?0—8)} and w_(A) = }{cotA@—,/(cot?0—8)}, (3.5) 
which are real if 9 < #, as we may suppose in this section. The paths of steepest 
descent are sketched in P, fig. 5, in the plane of a convenient variable 7 = sinh~1u 
in terms of which the function 
F(u) = (cos? —usin@) ./(1+u*) (3.6) 
is single-valued. The paths of steepest descent C_(4), C_(@) in the cut u-plane are 
their conformal images and are shown schematically in figure 1. 
The path of (3.4) can be deformed into C_+C.,. Consider the behaviour of 
the integral as 9 > 0. It is easy to see that the integrand tends to zero uniformly 


as |u| > 00 along C'_(9), and so, as # + 0, 
(l+w)exp{iN /(1+u?)} du 
(l+w)exp {iN /(1+u?)} du, 


whence im | -—-H{(3Y (Ax) — ¥,(A2)}, 


see equation (3.3)). 
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As for the integral along C_,.(#), we note that w,(@) ~ }cot@— oo; thus C,(A) 
moves off to infinity as 0 > 0. Introduce a new variable v = wtan 0, in terms of 
which 


* 0 exp(—tni) 


iN os? @ * tan- € ° 
(3. 7) 


The expression { } has a saddle point at v, (9) = w,(0) tan @ = 3(1+ ,/(1—8tan?@)), 
which remains finite as 9 > 0; and the expression (N cos? @/sin @) = M is a large 


= cot 0 (v2 + tan? ovexp| 


J ¢ x exp ($73) 








FicurE 1. Deformation of the contour of integration for small values of @. 


parameter tending to «0 as 0 + 0 (whether N is large or not). If we apply the 
ordinary method of steepest descents (see, for example, Jeffreys & Jeffreys 
1946, §17.04) to (3.7), we find (the details are omitted) that 


F cot? 7 


}. ~ m exp {iMf(#)—47i}, whence near 0 = 0 


acs Pkt Co(x.y) = im [ 


(NV cos? @ 


scene ey a Po 
= 172 N-2 (sin @)-? sin - 
; | sind 


f(0)— 3m +0(N-3(sin @)-%), (3.8) 


f(A) = g(t —},/(1—8tan?0))? (1+ ./(1— 8 tan? 4))}, 


>i as 0-0. 


The equation (3.8) shows that the free surface oscillates with infinitely increasing 
amplitude and infinitely decreasing wavelength as 7 > 0. This behaviour may be 
compared with the behaviour of waves due to a concentrated initial pressure 
(Lamb 1932, §239), and is responsible for the infinite wave resistance of the sur- 
face pressure point. (For a submerged pressure point or distributed pressure 
the amplitude on # = 0 and wave resistance are finite.) We note again that to 
obtain (3.8) it was not necessary to assume that N is large. 
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4. The neighbourhood of the critical lines 


We have noted in §1 that existing expansions are inadequate in various ways 
near 6 = +6,. On the line 7 = 0, itself, the function F(w) of (3.6) has a double 
zero, and the saddle points w,(4),u_(@) coalesce into a saddle point of higher 
order at u = 1/,/2 which is easily treated by the method of steepest descents. 
Peters (equations (5.3), (5.5) with an error in sign in (5.3)) has shown that 


7Pg 
——~—, C(x, x tan 0,) 
P, K* . 
3rd) . (34 .)\ 3 ‘oF 
~— “sin ( -N)- 22) cos ( “WJ, (4.1) 
23 N3 ZY 23 N32 Z 
where N = Kzxsec6.. 


Thus, on (and presumably near) @ = ,, the amplitude decreases as N-*, while 
elsewhere it decreases as N-* (see P (5.2), (5.4)), and so the waves are prominent 
near 7 = #, when_N is large. We shall now obtain an asymptotic expression which 
is valid in a finite angle inciuding the line 7 = 6,. This is to be compared with 
Hogner’s equations (H (82), (83)) which change in form as the line @ = @, is 
crossed; it is believed that these hold when N?|@—6,| is large and N?|0—0,| is 
small, but their precise region of validity has not been investigated. Our asymp- 
totic formula will be derived by the method of Chester et al. (1957) and will 
involve the Airy function (see Jeffreys & Jeffreys 1946, § 17.07) 


Ai(Z) = exp (4t3 — Zt) dt, (4.2) 


] foc exp (4z7i) 
= 

-710 . 2 exp ( ni) 
i 

= | cos (4t° + Zt) dt 

7 J 0 

when Z is real. Since in (4.2) the argument of the exponential is a cubic poly- 
nomial, comparison with (3.4) suggests the introduction of a new complex vari- 
able » by the implicit relation 


F(u, 0) = (cosO—usin@) (1+ u*) = — }e3 + (A) v—v(8), (4.3) 


where (4) and v(@) are to be chosen so as to make the transformation regular in a 
uniform neighbourhood of u = 1/,/2 for all 6 sufficiently close to 6,. The segments 
of the path of integration of (3.4) lying outside this neighbourhood are deformed 
into segments of paths of steepest descent from the saddle points, and the 
contribution from these is negligible, as in the ordinary method of steepest 
descents. The only significant contribution thus comes from the neighbourhood 
of the (nearly coincident) saddle points where the transformation (4.3) is 
applicable. 
On differentiating (4.3), keeping @ fixed, we get 


OF (u. é 
F(u 7) du oe. (4.4) 
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u dv 
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The left-hand side vanishes when uw = w., (9), the right-hand side when v = + 2(6). 

If the (w, v)-transformation is to be regular, these points must correspond, and 
so, from (4.3), 

F(u.,9) = (cosO—w, sin9)/(1+u?_) = 343(A)—v(4), ) 

F(u_,™) = (cosO—u_sin 6) /(1+u2) = —33(0)— v(9),) 

where w,(#) and u_(#) are known functions of 0, see equation (3-5). These equa- 


tions are readily solved for “(@) and v() 


33 cos? 0 ‘ 


180) = Gang (+ OF (1-30)? - (1 Q) (1 + 4Q)4, (4.6) 
3 2 ‘ 
(0) = — > 9 1 4 yh (1 —2Q)F+ (1 —Q)E( + 4Q)H, (4.7) 


32 sind ' 
where Q = ,/(1—8tan*@). By expanding these expressions it can be shown 
that (9) and v(@) involve only even powers of Q and are therefore regular func- 
tions of 0. 

This also follows from the general theory (Chester et al. 1957), which shows 
further that with these values of (#) and v(@) the transformation (4.3) is indeed 
uniformly regular, and that the points vu = w.(@),v = + (0) correspond. Follow- 
ing the procedure explained by Chester et al. (1957) we change the variable of 
integration in (3.4) to v, and expand (1 +?) du/dv in the form 


1 oO ve 
(1+ 08) T= SPn(8) (0? — (6) " + 2S. dnl) (02 — (8), 
) 0 


which holds uniformly when v and 4 — @, are sufficiently small. The theory shows 
that the asymptotic expansion of (3.4) is 


“pm(A) | (0? 1(9))™exp [iN ( — Jo + pw — v)] dv 


. 
a 


+ Xq,,(9) | u(v? — u())" exp [1N( — 403 + nv — v)] do, 
where the integration can be extended with a negligible error from —coexp (§771) 
to coexp(—4nz), ef. P, fig. 6. To obtain the dominant terms it is sufficient to 
calculate the leading coefficients p,(@) and q,(@). On putting w= wu, (#),v= + 13(0), 
we find that s ; 
(l+u2 )(<*) = po(9) + 139) q9(9), 


and the left-hand side is known when du/dv is known at v = + 04(0). 
But, from (4.4), 


oF a) OF d*u 
ae 


ou? \dv} © du dv? — 
(O2F ‘du\? 
whence ( =) (: = = 28, (4.8) 
cu?) . \dv, 


where (: 4 = sind QU. + Q)-*(1 F 1Q)-4, (4.9) 
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as is easily shown. The sign of (du/dv), is positive for small 6—6,, since u = u., 
vy = +? correspond. Thus, from (4.8) and (4.9), 

(du t cos? O / | (9) t . ‘ 

l+u2 = 1+Q)'(1510)8, 

' dv, Qe oe 1—8tan?0 (LEQ) (1+ sey, 


whence 


3! cos? / H(A) 
— aN _19)% oe oe 
P19) = a (i HO) 4 (1+ Q)# (1 QPS (1—Q)E (1+ 4Q)H 
(4.10) 
and ; ; 
34 cos?0/ (A) a 5 — i : 
a0?) = 2 sink @ | 1— Sor Q-M(1+Q)E (1 4Q) — (1 - QF 1 + 30) 


(4.11) 


where Q = N(1—8tan?@); the functions p,(9) and q,(9) are regular near # = 0. 
We thus see that the integral (3.4) is asymptotically 
pol) | exp [tN (— 403 + po rMde +460) vexp [tN (— 403+ wo—v)]dv 
= —ip,(A) } exp [V(iw3 + ww —1v)] dw —qy() | wexp [N (iw + ww —1y)] dw, 
where w = iv and the integration is from w = cexp(—4771) to w = oexp (477). 
But these integrals are multiples of the Airy function (4.2) and of its derivative. 
Thus the integral in (3.4) is asymptotically 


2mi exp (—iNv(A)) {- in-4 p(9) Ai(— Nin(O)) +N -49,(0) Ai’ (- Nin(6))}. 


Higher terms in the expression are of order N-! Ai, N-? Ai’ at most. Theamplitude 
C(x, y) is obtained from the imaginary part, ef. equation (3.4) 


PY +t» Polé /) _ ne , 
SP, K2S2"¥) ~ Nn} Ai (— N5(@)) sin (Nv(A)) 
0 i 
™ tai (—V%4(A)) cos (Nv()), (4.12) 


where (9), v(P), Po(P), Yo(P) are defined by equations (4.6), (4.7), (4.10), (4.11), 
respectively, and are regular near / = @,. It can be shown that near 0 = 6... 


3 
1(0) = —— (0—6,) + O((8 —6,)?), | 
93 
a a (4.13) 
v(8) = — = rel (0 —0,) +0 ((0—0,)2),| 
2 os 
33 15 
and that P,(9.) =, 99.) = = 
»3 93 


Values of the functions (4), »(9), po(P), Yo(A) are shown in Table 1. The functions 
Ai and Ai’ are tabulated by Miller (1946). 

Equation (4.12) gives the wave pattern far behind the travelling pressure 
point; from the general theory (Chester et al. 1957) it follows that it is valid in 
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0 (degrees) HG) v(@) P(A) Qo(P) 
10-00 0-532 — 1-227 9-484 11-131 
10-25 0-510 — 1-210 8-890 10-545 
10-50 0-488 — 1-193 8-348 10-002 
10-75 0-468 —1:177 7-852 9-500 
11-00 0-448 — 1-162 7:397 9-035 
11-25 0-429 — 1-147 6-979 8-603 
11-50 0-410 — 1-133 6-594 8-201 
11-75 0-392 — 1-120 6-239 7-826 
12-00 0-374 — 1-107 5-911 7-476 
12-25 0-357 — 1-095 5-670 7-148 
12-50 0-340 — 1-083 5-326 6-842 
12-75 0-324 — 1-071 5-065 6°555 
13-00 0-308 — 1-060 4-822 6-285 
13-25 0-293 — 1-050 4:596 6-032 
13-50 0-278 — 1-039 4-385 5-793 
13-75 0-263 — 1-030 4-189 5-569 
14-00 0-249 — 1-020 4-005 5:357 
14-25 0-235 —1-011 3°833 5-157 
14-50 0-221 — 1-002 3672 4-967 
14-75 0-208 — 0-993 3°521 4-788 
15-00 0-195 — 0-984 3°379 4-619 
15-25 0-182 — 0-976 3°246 4-458 
15-50 0-169 — 0-968 3-120 4-°306 
15-75 0-157 — 0-961 3-002 4-161 
16-00 0-145 — 0-953 2-891 4-023 
16-25 0-133 — 0-946 2-786 3°892 
16-50 0-122 — 0-939 2-686 3°768 
16°75 0-111 — 0-932 2-593 3-649 
17-00 0-100 — 0-925 2-504 3°536 
17°25 0-089 —0-918 2-420 3°428 
17-50 0-078 — 0-912 2-340 3°324 
17-75 0-067 — 0-906 2-265 3°226 
18-00 0-057 — 0-899 2-193 3-131 
18-25 0-047 — 0-894 2-125 3-041 
18-50 0-037 — 0-888 2-060 2-955 
18-75 0-027 — 0-882 1-999 2-872 
19-00 0-018 — 0-876 1-940 2-793 
19-25 0-008 — 0-871 1-884 2-716 
19-50 — 0-001 — 0-865 1-831 2-643 
19-75 — 0-010 — 0-860 1-780 2-573 
20-00 — 0-019 — 0-855 1-732 2-506 


TABLE 1. The coefficient functions “(9), v(@), 7o(7), Yo(7) occurring in equation (4.12) 


for the elevation 


some finite angle including 0 = and it reduces to the equations of Peters (P (5.2), 
(5.4)) when |N*(4)| is large, as can be shown by using the asymptotic expressions 
nas sana toe 2 awl 4 - 1 +>} ees 
Ai(X) ~ 37-2? X-* exp(—2X?), Ai’ (X) ~ —40-*X* exp (— 2X%),) 
; ‘ r oy} . , Ly 23 
Ai(— X) ~ 73X-4cos(2X!—4m), Ai’ (—X) ~ 1-?X! cos (2X? — 37), J 


3 4 


(4.14) 


valid when_X is large and positive. The equation (P (5.5))isobtained by supposing 
that N#(0) is small, expanding the Airy function in powers of @ — 0, and retaining 
only the first two terms. 
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5. Description of the wave pattern 


Let us consider first the well-known pattern some distance inside the critical 
lines where the Airy functions may be replaced by their asymptotic expansions 
in terms of circular functions. When _N is large, the surface elevation in this region 
is of the form (see P (5.2)) 


A,(9 me A,(O P 
cos {Nf,(0) +64} 4 - cos {Nf,(0) + ¢}. (5.1) 
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FiagurE 2. Contours of equal magnitude of the dimensionless surface elevation 
pU%E, (x, y)/2gP, due to a concentrated force P, moving with velocity U, for gr/U? = 2mr/A 


near 100. (Here r denotes the distance from the disturbance, A the wavelength of waves 
moving with phase velocity U.) For the sake of clarity only positive values of the eleva- 
tion are shown. Note that the largest elevation in the figure occurs near the point 
(120-2, 28), not near 6 = @.. 


The curves of constant phase corresponding to each term are of the form 

rf(@) = const.; 
that is, they are geometrically similar with respect to the origin (see also, Lamb 
1932, $256). The amplitude of each term falls off like r + The total surface 
elevation, being the sum of two terms, does not exhibit exact similarity with 
respect to the origin, but the amplitude is almost periodic along radii from the 
origin. 

Next let us consider the surface elevation (4.12) near the critical line 6 = @, 
when N is large. For very large NV and bounded values of Nilo- 9,| the first term 
in (4.12) is dominant. Thus, the nodal lines (contours of zero surface elevation) 
lie near those curves where the first term vanishes, that is, near those curves 
where either sin (N'(9)) or Ai(—N3,(0)) vanishes, that is, 


near Vv(?) = m7, where m is any large integer, (5.2) 


and near Nin(@) = isl, Ng. «25 Mile ~~<y (5.3) 


; 
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where a, < 0 is the sth zero of Ai(z). It follows from (4.13) that the curves (5.2) 
are ultimately equidistant straight lines crossing 0 = 0, at an angle tan-1 1/,/2 
that is, inclined at an angle tan-!,/2 = 37+ 40, = 54 3°, to the track of the dis- 
turbance, as was stated by Kelvin (W. Thomson 1891, p. 485). As for the curves 
(5.3), it follows from (4.13) that they are approximately of the form 
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FicurE 3. Contours of equal magnitude of the dimensionless surface elevation 
pU%,(x, y)/2gP,, for gr/U? near 1000. For the sake of clarity only positive values of the 
elevation are shown. The elevation again exceeds 0-08 near the point (963, 285). 
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where p is the perpendicular distance from 0 = 4, and r is the distance from the 
origin. These curves lie inside the critical lines, and it is seen that p increases only 
slowly with distance from the origin. Thus, the curves (5.2) and (5.3) combine to 


form a net of (approximate) parallelograms: the sides parallel to 0 = 0, are of 


nearly constant length, while the length of the sides parallel to 0 = 41+ 44, 

increases as the cube root of the distance. The amplitude due to the first term in 
: ° ° yi a ae . 

(4.12) varies with distance as V~* and this is ultimately much larger (by a factor 
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FicurE 4. Contours of equal magnitude of the dimensionless surface elevation 
pU%G, (x, y)/2gPo, for gr/U? near 10,000. For the sake of clarity only positive values of the 
elevation are shown. Note the progressive lengthening of the crests with distance from the 
disturbance, as shown by figures 2 to 4. 
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N') than the term N-! describing the variation with distance well inside the 
critical angle. Thus, ultimately the largest waves are near the critical lines 
= +0, and near 6 = 0 where the linearized theory predicts infinite amplitude 
and infinitesimal wave length. 

Figures 2, 3 and 4 show contours of equal surface elevation computed from 
(4.12) for three distances, near V = 100, 1000 and 10,000, respectively, where 
N = gr/U* = 2nr/A, and A is the wavelength of a regular wave-train travelling 
with phase velocity U. In each case the nodal lines lie near the boundaries of the 
parallelograms mentioned above. There are striking phase changes near the 
zeros of Ai(—N',(0)); these arise from the presence of the second term in (4.12). 
The lengthening of the crests with increasing N is clearly shown. Figure 3 
(N = 1000) resembles closely the contour map computed by Hogner (1923). 

In figure 2 (N = 100) the largest amplitudes occur near the point (120-2, 28) 
which is not near @ = 7, where they would be expected. The reason is that the 
quotient N- obtained in the last paragraph but one is about 0-46 and differs 
little from unity. The reduction in wave height by this factor is more than 
compensated by the fairly rapid increase in p,)(@) and q)() as 7 decreases from 
0,; see Table 1 above. Even when N is near 10,000, the quotient N-* is still about 
(0-22. Thus the concentration of amplitude near the critical lines is not very strik- 
ing (much less so than is suggested by Havelock’s comparison (1908, Table 2) of 
the amplitude on the critical lines with the amplitude of transverse waves on 
(= 0). 

In the computations it is assumed that the pressure is concentrated at a point 
on the free surface. If the pressure is distributed over an area, the resulting ampli- 
tude is obtained by an integration over the area, and in regions where the wave- 
length is short (compared to the dimensions of the pressure area) the amplitude 
will be reduced by destructive interference. In particular the infinite amplitude 
predicted near the track of the disturbance will not be observed. Similar remarks 
apply when the waves are due to a submerged instead of a surface source. 


Mr H. P. F. Swinnerton-Dyer computed for me the co-ordinates of points on 
the contours of equal elevation which are shown in figures 2, 3 and 4. His help 
and advice is most gratefully acknowledged. 
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Magnetohydrodynamic flow of a viscous fluid 
past a sphere 
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The flow of a viscous incompressible electrically conducting fluid past a sphere 
is studied; the uniform ambient flow field is colinear with the ambient uniform 
magnetic field. The force exerted on the sphere is computed for various con- 
ductivities and Reynolds numbers; of particular interest is the distinction in 
behaviour between the flow with ambient particle speed greater than ambient 
Alfvén speed and that with particle speed less than Alfvén speed. 


1. Introduction 

The flow of a viscous incompressible fluid past an obstacle at low Reynolds 
number has been the subject of many investigations. The rigorous analysis of this 
problem requires the solution of the non-linear Navier-Stokes differential equa- 
tions. Stokes (1945), Oseen (1910), and Lewis & Carrier (1949) have formulated 
linear problems whose relevant solutions are good approximations to the observed 
physical facts. The linearized problems have been solved for flow past a sphere by 
Lamb (1945) and by Goldstein (1929) and the drag on the sphere has also been 
found. We shall use the same methods of linearization to study the flow of a 
viscous, incompressible, and electrically conducting fluid past an obstacle when 
the magnetic field and velocity are constant and parallel far from the obstacle. 
A detailed solution will be obtained for a sphere and a method of solution for 
other geometries will be indicated. The modification introduced in a uniform flow 
by an externally applied point force will be described. 


2. The flow past a sphere 

The magnetohydrodynamic flow of an incompressible viscous electrically 
conducting fluid of constant properties is governed by Maxwell’s equations and 
the laws of conservation of mass and momentum. In m.k.s. units these equations 
take the form 


p . +IV.v)V4 = —Vp'-+ py¥2V' + yj x H’, (2.1) 
a 

div V’ = 0, . (2.2) 

curlH’ =j, divH’ = 0, (2.3) 

curl E = 0. (2.4) 


These are supplemented by the constitutive relation 


j = o[E+yzV’ x H’]. (2.5) 
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Here, V’ is the fluid velocity, H’ is the magnetic field, E is the electric field, is 
the electrical conductivity, p is the mass density, is the magnetic permeability, 
and pv is the kinematic viscosity. The differentiations are taken with respect to 
the physical variables 2’, y’, z’. 

The differential equations can be put in a dimensionless form by introducing 
the following substitutions: 

a cd y’ 4 H’ 
V= TW’ ro ae te, Pe, te and H= Hi 

Here, U is a characteristic velocity of the problem, H) is a characteristic magnetic 
field, a is a characteristic length and the differentiations are with respect to 
x, y and z. In the problems we are considering U is the undisturbed velocity 
and H, is the magnetic field far from the obstacle. For steady-state problems, 
éV/ot = 0 and (2.1) to (2.5) take the form 


2 
R(V.V)V = = Vp + ev x H, (2.6) 
div V = 0, (2.7) 
curlH = ae div H = 0, (2.8) 
Ay 
curl E = 0, (2.9) 
j = o{E+yUH,V x H}. (2.10) 


If we substitute (2.8) into (2.6) and use the vector identity 


(curl H) x H = (H.V) H—3V(H.H), 


we obtain 


M? 
R(V.V)V = -VpytVV+5 (H.V)H, (2.11) 
M? : 
where Do = pt oR H.H. (2.12) 


Here, R,,, is the magnetic Reynolds number given by R,, = Uayo and M is the 
Hartmann number given by M = ~H,a(c/pv)?. The ratio R,,/R is a measure of 
the relative importance of magnetic and viscous effects. 

If we take the curl of (2.10) and make use of (2.7) and (2.9) and the vector 


meet curl (curl H) = V(div H)— V?2H = —V?H (2.134) 
and curl (V x H) = (H.V) V—(V.V)H+V(V.H)—H(V.V), (2.138) 
we obtain V7H = — R#,,{(H.V) V—(V.V) H}. (2.14) 


We shall confine our attention to problems in which the fluid flows past an 
axially symmetric obstacle in an otherwise unbounded domain with an ambient 
velocity i¥, and magnetic field iH, which are uniform and directed along the 
symmetry axis of the object. Chester (1957) has treated this problem following 
a Stokes-like analysis; he computed the drag on a sphere to be 


D = D{1+3M + sig M?—7435M*+0(M4)}, (2.15) 


where D. = 6npvaU. (2.16) 
28 Fluid Mech. 8 
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We shall invoke an analysis similar to that of Oseen, although the motivation 
stems from the Carrier—Lewis approximation technique as used by Greenspan & 
Carrier (1959). In this technique one anticipates that the mathematical model 
retains significance when each of the undifferentiated quantities in the convective 
terms of equations (2.11) and (2.14) is replaced by an appropriate average which 
is taken to be some fraction of the free-stream values. In this paper, we shall 
take these functions to be unity, so that we shall treat a mathematical problem 
which is identical with that which would arise if the Oseen philosophy were 
followed. It should be noted, however, that the equations which would be 
obtained with a more appropriate choice of these averages would have solutions 
which are related to those obtained here by the formulas 


Y Y R R,, C - 
V,(r,C,, C,, C,, 8, R,,, 8) = V(r. 24, Lo C,” af), (2.17a) 

an R R, C , 
Hy(6, Cy, Cp,Cy, By RB) = ACE LL AG oe Ze), (2.176) 


where Ci, C,i, C,i, C,i are the quantities which replace V, H, H, V, respectively, 
as they appear in (2.11) and (2.14). 
With the foregoing, (2.11) and (2.14) become 
OV M?cH 
R—=-V vey ~— 
Cx Pot ig R,, 02 


m 


V?H = —(H-V). (2.19) 
Cx 


(2.18) 


Since the problem has been linearized superposition may be used; we write 
V =i+vand H =i+h. The boundary conditions require that v and h tend to 
zero as the distance from the obstacle goes to infinity, that v = —iat the obstacle 
surface, and that A, and wh, be continuous at the obstacle surface. Here, h,, and 
h, are, respectively, the component of magnetic field that is normal to the 
obstacle and the component that is tangent to the obstacle. 

Equations (2.18) and (2.19) can be simplified in the following way. We define 


Cy ; 
Vv, = v——-h, 2.20 
: Rk, ( 
ile 
v, = v——h, 2.21 
a salt 
where the «; are the roots of 
a’?+(R-R,,)a—-M?=0 (2.22) 
L . CV 
We then obtain R, a = —Vp,+ V2v,, (2.23) 
a 
CV, 
R, 5 = —Vp,t+V2Vo, (2.24) 


where Rk; = R+a;. Each of these equations is identical with that encountered 
in the classical Oseen treatment of viscous fluid problems and our investigation 
will parallel closely that of Oseen. 
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It is easily seen that 


] 
R 
= —— fy, 2.2 
h R,-R, [Vy Vol, ( 6) 
div v, = divv, = 0, (2.27) 


and that the ‘modified’ Reynolds numbers R, are the roots of 


R?—-(R+R,,) R;+(RR,, — UM?) = 0, (2.28) 
and are always real. 
2 2 
We define Poe Pe (2.29) 


pU2 RR 


m 


and note that, for # < 1, each of the ‘modified’ Reynolds numbers is positive, 
but for # > 1 there is a positive and a negative ‘modified’ Reynolds number. 
These two cases will be considered separately. When # = 1 the ‘modified’ 
Reynolds numbers are R + &,,, and zero. 
It is convenient to follow Lamb (1945) and introduce the potential ¥; such that 
, ee ™ : on(r-*) 
V, = —X;i+ Rj gradX;+grad > A, - ; (2.30) 


ae 
n=0 J 02 


In order that (2.23) and (2.24) be satisfied it is necessary that 


|"- R, =| xX; = 0. (2.31) 
00 Ant+1(y—-1 
Furthermore, Po=— > A, R;* : P na ) (2.32) 
n=0 . 


The axially symmetric solution to (2.31) that vanishes at infinity is 


X; = 2 Y Bi Xa(|Ri| 7) P,(cos 4), (2.33) 


n=0 


where @ is the polar angle measured in spherical co-ordinates (r, 0, w) and 


} 
x,(x) = (2n+1) (") Kya(5): (2.34) 


Here, K,, is the modified Bessel function of the second kind. 

It follows from (2.8), (2.26), (2.30) and (2.31) that the current paths are circles 
which have their centres on the z-axis. From the axial symmetry of the problem, 
equation (2.10), and Kirchhoff’s law, we deduce that E = 0 everywhere in the 
fluid. Inside the spherical obstacle the electric and magnetic fields are governed 
by Maxwell’s equations and the charge conservation equations which, in m.k.s. 
units, are given by 


curl h; = divh,=0, j,= o;E,,| 


ae 
pi 
curlE;=0, divo,E, = 0. 
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Since the electric field vanishes in the fluid and on the boundary of the sphere, 
it follows that E; = 0 and that 


h. = _—e, . D,,1r"P,,(cos 4). (2.36) 


L 5 
n = 


We now follow the method of Goldstein (1929) and write the boundary con- 
ditions at the surface of the sphere (r = 1) in the form 


(n+1)(—1)"R,, l = Rd ; 
n! Ro “wee. MY Bm Bp) Bim Fn 
(0 (sn = 0,2,3,..:)) joan 
- 2.3 
ll (n=1), , eee 
oe 1)" R,,, . l = = j 
mR ont oT eer ee lia (Bn — Bs) Bim Tm) 
_f09 (n=2,3 )) 
~~ \l (@w=1) a 
(—1)"(n+1)R R > = 
me / = m a B : 
n! R A” — 3(R,— Ry) nose PB im Fam) 
— TMD, Sn (n = 0,1,2,...), (2.39) 
l 
(—1)"R R - 
™A,+ ” —l) R,)=D,+6 ee 
n! R Ant o(R,— Bi) nwo | YB inTn,m(H5) n tin (a 
(2.40) 
where 
on, On mR, r) = 2{v,,(|Ry| r)y nim (\R, |r) age Xm (|R,| 7) Wn, | (|R,| r)} (5 )’ 
Ry 
(2.41) 
4 = : 2X m(|R,| ") 2X By »? 49 
and Tn, m( yr) 7 R, hy mi | R,|) = n,m (|R,| 7) | R ‘ (2. <) 
The functions X,, ,, and ‘¥’, ,, are given by 
‘ _ m(2n +1) (2m)! (2n)! — 
“nm n+l (2n+2m)! \(n)!(m)! 
(Wirim-1, 3(n—1) m—1 2n+2m—-1 Vntm—3 , 1.3.5 (n—1)(n—2) 
“| £ “(Qn—1)2m—-1 n+m E 2!) (2n—1)(2n—-3) 
— (m=—1)(m—2) (2n+2m—1)(2n+ 2m—3) Prim-s 
* (2m —1)(2m—3) (n+m)(n+m-—1) ‘2 
a 1.3... (2r+1) (n—1)...(n—r) (m—1)...(m—r) 
ere r! (2n—1)... (2n—2r4+1)(2m—1)...(2m—2r+1) 


; (2n+2m—1)...(2n+2m+4+ 2n4 1) Wysm—or-1 
(n+m)...(n+m—r-+l1) 


eo (n—1)...(n—m41) (2n+2m—1)... (2n+3) Pr_miil (2.43) 
(2n—1)...(2n—2m+3) (n+m)...(n+2) g& | 
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and 
Y= (2m)! (2n+1)! /(n+m)! ) | , MN m 2n+2m+1 y 
nm “~ (2n+2m+1)! lane \ nem In —12m—1 n+m aes 
12.3 n(n — 1) m(m— 1) (2n + 2m+ 1) (2n+ 2m— l) 
ay (2n — 1) (2n —3) (2m — 1) (2m— 3) (n+m)(n+m—1) ¥n+m—4 
_ 1.3.5... (47-1) n(n—1)...(n—r+1) 
aor r! (2n—1)(2n—3)... (2n—2r+1) 
: m(m—1)...(m—r+1) 
* (2m — 1) (2m—3)... (2m —2r+1) 
(2n+2m+1)(2n+ 2m—1)... (22+ 2m—2r+ 3) y 
Sn+m—2r 


(n+m)(n+m—1)...(n+m—r+l]1) 


. n(n—1)...(n—m+1) 
. “*T (In —1)(2n—3)... (2n—2m4 1) 


(2n+2m+1)(2n+ 2m—1)...(2n4+3) , 


x 2.44 
(n+m)(n+m—1)...n+1 Wn > ( ) 


: m\% x , 
where r,(") = (2n+1) 4 | i( ). (2.45) 
Equation (2.37) is the condition v, = — cos @, equation (2.38) is the condition 


vy = sin#, equation (2.39) expresses the continuity of wh, and equation (2.40) 
represents the continuity of h,. The infinite system of equations is approximated 


to order q by considering (2.37) for n = 0,1,...,q; (2.38) for n = 1,..., q; (2.39) 
for n = 0,1,...,¢g—1; and (2.40) for n = 1,2,...,q. This yields 4q¢+ 1 equations 
for the 4¢+1 unknowns Ao, Aj,...,4,; By, By, ---,Big13 Boo, Bay, ---, Bogs; 
D,, D,, ..., D, with the rest of the A,, B;,_;, D, taken as zero. Numerical results 


for the constants were obtained using a second-order approximation and 
several values of the ‘modified’ Reynolds number. For small ‘modified’ Rey- 
nolds numbers it was found that Ay, A,, By, By) and D, are all much greater than 
A,, B,,, By, and D,. 

The drag D on a solid cbstacle is obtained by applying Newton’s second law 
of motion to the fluid surrounding the obstacle and bounded by a surface S 
which is everywhere far from the obstacle. The total force exerted on the fluid 
within S is the sum of the viscous force exerted by the obstacle (—D), the force 
exerted by the fluid outside S, plus the ponderomotive force. 

We are able to show by direct calculation from the results of § 2 that for small 
‘modified’ Reynolds numbers the ponderomotive force {#R(curl H) x H} is 
small compared with the viscous forces vVV in all parts of the fluid provided that 
fu = p,. This restriction on the permeabilities is necessary so that near the sphere 
H ~z i. Hence in computing the drag we neglect the body force. This calculation 
can be carried out by an extension of the process used by Lamb in showing for 
small Reynolds number that (curl V) x V may be neglected compared with viscous 
forces in a non-conducting fluid. 
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We continue in a manner similar to that used by Goldstein (1929) and obtain 


D R 0d 
Cp = sas =~  (I- ~-d8. 2.46 
D”™ pU?a R | p)| On aid 
For flow past a spherical object we obtain 
4nR 
Cy = EM(1-f) Ay (H = M)- (2.47) 


The first-order solution for small ‘modified’ Reynolds numbers yields 


677 




















Cp & RP (1+2R4+O0(R?)) (P< 1), (2.48a) 
67[, 3R/R,,-2R,,fB-R\ 
dQ ~ a m m > : 2. R 
when s/t = /,. 
16 T T T T 
5b 
l4b 
13 
§ 
E iak 
Oe 
10F 
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FicureE 1. Drag coefficient versus 2. ¢€= R,,/R, R = 0-5. 


Numerical values were obtained for C, at R = } for several values of R,, and 
£, where all terms up to those of order R? were retained and it was found that 
(2.48a) and (2.485) were changed by about 1°%. These results are plotted in 
figure 1. This does not mean that the error in these formulas for Cp is as small as 


1%, since this error is dependent on the accuracy of the linearization. 

In order to obtain additional insight, the results of §2 were used to obtain 
numerical values for the x-components of particle velocity and magnetic field 
along the lines (0,0,z), (5,0,z) and, for the subsonic case, (—5,0,z), with 
R= R,, = 0-25 and wv = ,. The results are shown in figure 2. It was found that 


both the particle velocity and magnetic field are independent of # for 0 < # < 1. 
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In fact, in both the supersonic and the subsonic case (f = 2), the magnetic field 
was found to be small (less than 0-1) and hence was not plotted. If the sphere has 
a permeability different from that of the fluid, the magnetic field is distorted 
near the sphere but the effects are unimportant elsewhere. 

As was mentioned previously, Chester (1957) extended the Stokes-type 
approximation for the flow past a sphere. This is equivalent to the linearization 
we have used, if the special case R and R,, > 0, # large and yw = p, is considered. 
This yields ‘modified’ Reynolds numbers of + M and a drag, particle velocity 
and magnetic field which are all in agreement with Chester’s results. 
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FIGURE 2. x-component of particle velocity versus position; R = R,, = 0-25. The positions 
of the origins of the three graphs indicate the values of x at which u(x, y) is plotted as a 
function of y. —, 8B < 1; ---, B = 2. 


3. The point-force problem 

We shall discuss here the flow field and magnetic field which arise when an 
externally applied point force acts on a fluid which otherwise would have 
uniform velocity and magnetic field. Such a description is of interest because it 
will be a major contribution to the description of the flow past objects such as 
that considered in § 2, and because it can be displayed in closed form with easy 
interpretation. 

The governing equations are (2.14) and (2.11) modified only by the inclusion 
of the point-force contribution. That is, 


M? 
R(V.V)V+Vp,—-V?V + = (H.V)H = — Fd(r)i, (3.1) 
with d(r) = d(x—0, y—0,2—90), 


where 6(r) is the Dirac delta function and — F'd(r)i is the force per unit volume 
applied at the origin. The equations of motion are put into dimensionless form 
and linearized, and linear combinations are taken in the same manner as in § 2. 
In this problem we define the characteristic length by 


(6zpvU 
a= [eo ) (3.2) 
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The exact solution to the linear problem is found to be 


y — 3exp[-$ Ril r+ MRD] 3 v( mela bale + HBA) 3 v(!), 
j J 


J 2 r oR. Tr / 2R, ie 
(3.3) 
30 /1 
and Po=t5 a (*) ‘ (3.4) 


Expressions for v and h are obtained from (2.25) and (2.26). We note that at the 
origin the velocity becomes infinite as r~! and is directed along the positive 
x-axis, while h = 0. 

At this point it is appropriate to compare our results for # < 1 with those for 
B > 1. The wake for the magnetohydrodynamic flow is the region in which 
exp[—$|R,| r+ 4(R;) 2] is not small. In the wake the particle velocity and mag- 
netic field vary as r~! far from the origin, while in the region outside the wake 
they vary as r-*. For 0 < # < 1, the wake appears near the positive z-axis. 
However, for # > 1 the wake appears not only near the positive x-axis, but also 
near the negative x-axis. 

We now find it convenient to turn our attention to the co-ordinate system in 
which the fluid is at rest at infinity. From (2.25), (2.26) and (3.3), it follows that 
far from the origin at points outside the wake for # < 1 the velocity and magnetic 
field are as if due to a source located at the origin and of strength 67R#,,/R, R,. 
There is an equal inward flux along the wake near the z-axis. If # > 1, the 
magnetic field outside the wake is the same as that due to a sink at the origin, 
which is o‘iset by concentrated outward flows near both the positive and the 
negative x-axes; the velocity is that which would be due to a sink equivalent to 
the magnetic field sink and, in addition, there is an inward flux near the positive 
x-axis, both of which are balanced by an outward flux along the negative x-axis. 
The results for # = 1 are obtained by considering what happens in the limit as 
£8 > 1—- or 1+. These limits yield the same conclusions. There is a net influx of 
particles near the positive x-axis which is balanced by a source at the origin. 
The essentially different behaviour for 8 < 1 and f# > 1 may be explained by 


noting that, if # = 1, My 
in /6) Hp, (3.5) 


which is equal to the Alfvén wave speed. Hence, for # < 1 the free-stream velocity 
is greater than the Alfvén velocity and for # > 1 the free-stream velocity is less 
than the Alfvén speed. 


4. The infinite conductivity case 
We shall now discuss the flow past an obstacle of a fluid which has infinite 

conductivity. For 7 > 00 (R,, + 0) a possible solution to (2.7), (2.11) and (2.14) 

is h = v, which reduces the problem to the form 

é 


vera nif sly = grad pp, (4.1) 


V.v=0. (4.2) 

















mb 
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These equations are identical with the equations of motion for a non-conducting 
fluid flow with Reynolds number R(1—/) and pressure py. For £ > 1, we deduce 
that the disturbance velocity is equivalent to the disturbance velocity due to the 
flow of a non-conducting fluid which, far from the obstacle, is in the direction 
of the negative x-axis and which has a Reynolds number R(f—1). Hence, if 
we can satisfy the boundary conditions on the velocity in the problem in which 
there is no magnetic field, then we can satisfy the boundary conditions on velocity 
in the magnetohydrodynamical problem. The magnetic field satisfies the same 
boundary conditions at the obstacle (H = 0) as the velocity as long as the obstacle 
has the same permeability as the fluid. The magnetic field inside the sphere is 
zero and the electric field is everywhere equal to zero by the same argument as 
in the previous section. For # = 1, (4.9) reduces to the equation obtained by 
Stokes for a non-conducting fluid. Naturally the suitability of this model can 
only be determined by the degree of agreement with experimental results. 


The work reported in this paper is an extract from a thesis submitted in 
April 1959 in partial fulfilment of the requirements for the degree of Ph.D. in 
applied mathematics at Harvard University. 
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Gravity waves on the surface of an inviscid incompressible fluid of finite depth 
are considered. The waves are assumed to be periodic in time and in the hori- 
zontal direction. The surface profile, potential function, pressure and frequency 
of the motion are determined (to third order) as series in powers of the amplitude 
divided by the wavelength. It is found that the frequency increases with ampli- 
tude for depths less than a certain multiple of the wavelength and decreases 
with increasing amplitude for greater depths. Graphs of the surface profile and 
of the pressure as a function of depth are included. 


1. Introduction 

Gravity waves on the surface of a liquid are governed by non-linear equations. 
In the classical theory of such waves, the equations are linearized. Thus, the 
results of that theory vield the linear terms in the expansion of the wave motion 
in powers of the amplitude. We shall determine two additional terms in this 
expansion for standing waves in a liquid of uniform finite depth. Our results 
will then exhibit various deviations from those of the linear theory. For example, 
the period will depend upon the amplitude, the motion will not be sinusoidal in 
time nor in space, and the maximum elevation will not equal the minimum 
depression. These effects become more pronounced the larger the amplitude. 
They are described and depicted graphically in § 4. 

Similar results have been obtained for liquid of infinite depth by Penney & 
Price (1952). Our results for the wave profile and the period coincide with theirs 
when we permit the depth to become infinite. However, our analysis indicates 
that in the expression for the fluid pressure there are additional terms, indepen- 
dent of the vertical distance, which they have not taken into account. Airy, 
Stokes, Rayleigh and others have obtained analogous results for progressing 
waves. 

We formulate the problem in the next section and solve it by a perturbation 
method in § 3. Then we discuss the results in § 4 where some graphs are also given. 
In the Appendix we show how the method of Penney & Price can be modified 
to yield the same results, although it requires more labour than our method. 


2. Formulation 


Let us consider the time-periodic irrotational two-dimensional motion of an 
inviscid incompressible fluid bounded below by a rigid horizontal bottom and 


above by a free surface. We suppose the motion to be periodic in the horizontal 
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direction and symmetric about the vertical plane x = 0. Then we may confine 
our attention to the fluid between that plane and a parallel plane one half wave- 
length from it. Let A denote the wavelength, k = 27/A the propagation constant, 
kh the mean depth of the liquid, k~!z and k-'y the distances along the hori- 
zontal and vertical axes, respectively, g the acceleration of gravity, (kg)? w the 
angular frequency and (kg)-+w—t the time. In addition let a be a measure of 
the wave amplitude, the precise meaning of which will become clear later. Then 
we define ¢ = ka and let ek~1y(x, t) denote the elevation of the surface above the 
mean level given by the plane y = 0. Finally, we introduce the potential function 
eg'k-id(a, y, t). 

In terms of the dimensionless quantities which we have just introduced, the 
equations which govern the motion are 


Ag=0 in O<a<a, —h<y<en(z,t), (1) 
n+ oG,+ t6(92+9;) =90 on y=en(z,t), (2) 
Py =omt+ep Nf, on y= en(z,t), (3) 
cg 
~-=0 =0,z=27,y=—-h, 4 
an on 2 xL=7,Y 2, (4) 
| n(x, t)dx = 0, (5) 
J0 
Vo(x, y,t +27) = Vd(z,y,t), (6) 
| | : n(x, t) cost cosadtdx = 0, (7) 
JoJo 
ro Cn (2n 
| | d(x, y,t) cost cosx dtdxdy = 4n*(tanhh)}. (8) 
J-hJIJ0 





Equation (1) expresses the incompressibility of the fluid; (2) is the condition that 
the pressure at the surface, as obtained from Bernoulli’s equation, be constant; 
(3) is the condition that a particle on the surface remain on the surface; (4) asserts 
that the bottom y = —/h is rigid and that x = 0 and 2 = mare planes of symmetry 
of the motion; (5) is the condition that the mean surface is y = 0; (6) asserts the 
periodicity of the velocity components; (7) and (8) fix the phase and amplitude 
of the motion. Equation (8) shows that in the dimensional units, a is propor- 
tional to the Fourier coefficient of 9(2,t) with respect to the linearized surface 
profile sint cosx. In fact a is just the amplitude of the linearized surface wave 
motion. 
The pressure p(x, y,t) is given by Bernoulli's equation 
(p—Ppo) = —y— eu, — 3€7(62 + 92). (9) 
Pg ; 
Here p, denotes the pressure of the atmosphere above the fluid and p denotes the 
constant density of the fluid. 
The problem we consider is that of determining the dimensionless surface 
profile (2, t), the dimensionless potential function (2, y, t) and the dimensionless 
angular frequency w satisfying equations (1) to (8). We note that the solution 
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also depends upon the dimensionless depth h, which is proportional to the actual 
depth divided by the wavelength, and the dimensionless constant ¢€ = ka, which 
is proportional to the linearized amplitude divided by the wavelength. We shall 
solve the problem by determining the first three terms in the expansion of the 
solution in powers of ¢. Additional terms can be found by continuing our 
procedure. 

We shall find that, for certain values of h, the problem we have formulated 
does not have a unique solution. This lack of uniqueness does not refer to the 
arbitrary constant which may be added to ¢. These values of h are those for 
which the linear theory yields, for some harmonic (in space), a frequency (in time) 
which is an integral multiple of the fundamental frequency. According to the 
linear theory the frequency of the mth spatial harmonic is (n tanh nh). Thus 
we can make the problem have a unique solution by imposing the condition 


n tanh nh 


integer)? =o ees 0 
ak + (integer)" (n 3 ) (10) 


Of course the solution we shall obtain will be a solution for any value of h, but 
it will not be the only solution satisfying (1) to (8) unless / satisfies (10). 


3. Perturbation solution 
We assume that 9, ¢ and w have limits 7°, ¢° and w, as € tends to zero. If we 
set € = 0 in equations (1) to (8), we find that all but (2) and (3) are unchanged in 
form while these two become 
n° +0,9) =90 on y=), (2°) 


oy —Uy =9 on y=. (3°) 


The zero-order problem, with (2°) and (3°) instead of (2) and (3), is just the 
classical linear problem. It has the unique solution 


n° = sint cos2, (11) 
g = =< j cost cosx cosh (y+h), (12) 
w2 = tanhh. (13) 


We now assume that 7, ¢ and w have derivatives with respect to € at € = 0 and 
we denote them by 7, ¢! and ,. Then we differentiate equations (1) to (8) with 
respect to € and let ¢ tend to zero. In differentiating (2) and (3) we utilize the 
relation d ‘ 5 + 

de (x, en, t,e) = ; +(y+€n,) <| Q. (14) 


a 


Upon differentiating (2) and (3), but before setting e = 0, we obtain 
Ne tO, oP T ol Dr. si (7 + €7), ) Pty] a a 3(P2 + ??) am E919 re + (7 * €N,) Pry] 
—€¢[Pye+ (9 +E.) Pyy] on y=en(z,t), (15) 


Pye a (9 = €N,) Pyy — ON — ON, = Ve ?, + EN x. ?, tr: EN AP re - (y + €N,) Pry] 
on y=en(z,t). (16) 








an 


an 





CO ed ~-— 





Standing surface waves of finite amplitude 445 


When we set € = 0 in the differentiated form of (1) and (4) to (7), these equations 
retain their form with the superscript ‘one’ affixed to 7 and ¢. From (15), (16) 
and (8), we obtain 


+ Moh = —3(Pz)? + (Py)7]— 09-9? on y=, (2*) 
oy — oon = O.92—9° Po, +0,92 on y=, (3!) 
| | ‘ n(x, t)sint cosadtdx = 0. (8!) 

0/0 


To solve the differentiated equations for y!, ¢! and w,, we first insert into (21) 
and (3!) the zero-order quantities given by (11) to (13). Upon doing this and 
simplifying the results, we obtain 


m+ Wo Pt = 4[(w2 — Wo *) + (W2 + Wo 2) cos 2x — (3«2 + Wo ?) cos 2t 
a ¢ aed Wy . ee 
— (3w9 — Wo “) cos 2t cos 2a]+—sintcosx on y=0, (17) 
Wo 


| a 

Py —VN =— 5, Sin 2t cos2x+,costcosr on y=0. (18) 
“Wo 

Next we differentiate (17) with respect to t and eliminate 9} from (17) and (18). 

This elimination yields 


$1 + w2dh = }(303 + wo *) sin 2t + 3(W3 — wo *) sin 2t cos 2x 
+2w,costcosx on y=0. (19) 


To determine ¢! we expand it in a Fourier cosine series since by (4), 61 = Oat 
x = 0and x =7. In order that ¢! satisfy (1) and (4), we find that it must be of 
the form E 


P(x, y,t) = ¥& A,(t)cosnx cosh n(y +h). (20) 


howd 
n=0 


Insertion of (20) into (19) yields 


0} Ao, = } (303 + wp *) sin 2, (21) 
w5coshh.A,,+sinhh.A, = 20, cost, (22) 

w§ cosh 2h. A,,+2sinh 2h. A, = }(w}— wp *) sin 2t, (23) 
wgZcoshnh.A,,+nsinhnh.A,=0 (n= 3,4,...). (24) 


From (6) and (20) we see that all A,, with > 1 must be periodic in ¢ with period 
27. Then from (10) and (24) it follows that A, = 0 for > 3. From (10) and (23), 
we have 3 


ea (Wy — Wp *) Sin 2¢. (25) 


in 


The periodicity of A, requires w, = 0, and (22) yields 
A, = a, sint+/, cost. (26) 
Here x, and /, are constants so far undetermined. Finally (21) yields 
Ay = — qg(3wo + Wo 3) sin 2t + at + Bo. (27) 


Here a, and f, are also undetermined constants. 
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We now insert the above results into (20), use the resulting expression for ¢! 
in (17), and solve (17) for 7!. This yields 
ni = 4(w2 — Wo *) — Wp %o — Wo(%, Cost — 8, sin t) cos x cosh h 
+4 (w2 + wp *) cos 2x + 4(wo 2 — 3w9 ®) cos 2t cos 2x. (28) 
By applying (7) to 7! we find a, = 0. Application of (81) to 9} yields £, = 0. 


Utilization of (5) leads to 
) Oy = § (Wy — Wo *). (29) 


Now that all constants except the inconsequential £, have been determined, the 
solution »!, d! and w, is completely determined. It is given by 


n) = 1[(w2 — wp?) + (wp 2 — 304 ®) cos 2t] cos 2x, (30) 
: 3 
d! = By t+h(Wy— wo 3) t — 71g (39 + Wo °) sin 2t — Ran ae (Wy — @ *) 
x sin 2¢ cos 2x cosh 2(y+h), (31) 
w, = 0. (32) 


It is interesting to observe that 7! contains a non-constant term which is in- 
dependent of t and that ¢! contains t-dependent terms independent of x and y, 
one of which is not periodic in t. 

Let us now assume that 7, ¢ and w have second derivatives with respect to € 
at ¢ = 0, and let us denote them by 7”, 6? and w,. To determine them we differ- 
entiate equations (1) to (8) twice with respect to e and set € = 0. Equations (1) 
and (4) to (7) remain unchanged in form while (8) takes the form of (8), with the 
superscript ‘two’ affixed to 7 and ¢. The twice differentiated form of (2) and (3) 
is most easily found by differentiating (15) and (16). To solve the resulting 
equations we proceed exactly as before. Upon eliminating 7? from (2) and (3°), 
we obtain 


9 


22 + w3 Gi = %, cost cos x + a3 cost COS 3x + 5) COS 3f COS + X33 COS 3t COS 3x. 


; 33 
Here the constants «;; are given by - 
241 = 2Wy + 7g (— Dwg * + 129% + 3er9 + 205),) 
O13 = zg (Wo * — 5w_ + 2098), ; 
(34) 


Xa, = zg (— 909 ’— 62w5 3+ 31wp), 
As, = 733(9Wo ’ — 22Wp + 13). ) 
In solving for ¢? as before we find that a,,, the coefficient of cost cosz, must 
vanish. This yields for w, the result 
Wy = gs(99 7 — 12W9 3 — 39 — 209). (35) 


Then calculations like those previously given lead to the following solutions for 
ny? and ¢ 


ny? = b,, sint cosx+b,,sint cos 3x + by, sin 3t cos x + bg, sin 3t cos3x, (36) 
¢? = £,+/,3cost cos 3x cosh 3(y +h) + £3, cos 3t cos x cosh (y +h) 
+ $33 cos 3t cos 8x cosh 3(y+h). (37) 








rey 
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Here f, is arbitrary and the other constants are given by 








bi, = g3(3w5 + 6wo *— 5 + 29), ) 
bis = zoa(9Wo 8+ 27wo 4— 15+ w§ + 208), (38) 
bs: = reg(3wo * + 189 45), | 
bsg = z35(— 9p 127+ 3w9 8 — 3w9 4+ 1), 
] 
aia Ses 1 4 3 —9  #,.—-1 9w3 
and Ais iihansk 3h | + 3w9) (3w9 ° — 5w9 * + 209), 
Ba = 128 cosh 7, (20 * + 62w9 °— 31wp *), , (39) 
1 
Bs = japcaee (1+ 3w§) (— 9wo 13 + 22weq ° — 136 5). 
4, Discussion 
The finite Taylor expansions of €7, e¢ and w are 
en = €n°(x, t) + €2y1(a, t) + de3y?(x, t) + O(e?), (40) 
ef = 69 (x, y,t) +e2D'(x, y, t) + Je°S*(ex, y, t) + O(C4), (41) 
@ = wy + fe*w, + O(e?). (42) 


The zero-order solution is given by equations (11) to (13), the first-order solution 
by (30) to (32) and the second-order solution by (35) to (37). 
On the basis of (42) we see that the frequency w depends upon the ‘amplitude’ 
e of the wave motion, as is usually the case for non-linear systems. By examining 
(35) we find that w, = 0 at w, = 0-89 which corresponds to h = h* = 1-07. For 
h > h*, w. < 0 while for h < h*, w, > 0. Thus, for depths greater than h*, the 
fluid behaves like a soft spring, its free vibration frequency decreasing with 
increasing amplitude. For depths shallower than h*, the fluid behaves like a 
hard spring, its free vibration frequency increasing with increasing amplitude. 
From (40) and the expressions for 7°, 7! and 7”, we see that the surface is never 
flat. It is most nearly so when t = nz where n is any integer. Then 
en(x, nm) = 1E?(w2 + 2w9 2 — 36 ®) cos 2x. (43) 
The velocity vanishes throughout the fluid at ¢ = (n+ 4)m. At these times each 
part of the surface is either at its highest or lowest position. In particular, when 
n is an even integer, the crest is at x = 0 and the profile is given by 
én = [ - Ld (9wp §+ 6w9 *— 154+ su) cos x + e?(wo ? + 3w6 ®) cos 2x 


251 
+ 53,63(9w9 * + Gu ® + 30wg * — 16 + wh + 20§) cos3x. (44) 
The greatest rise occurs at x = 0, and it is equal to 


x 


xe (275 1? + 27wp § + 96we ® — 63 + 11w4 + 6w8). 


(45) 


ENmax = E+ fE2(Wo * + Bw °) + 
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When ¢ = (n+4)7, with n an odd integer, the crest is at 2 = 7. The wave profile 
at this time may be obtained either by reflecting that given by (44) about x = 47, 
or replacing € by —e in (44). A graph of the surface profile for ¢ = 0-05 and 
h = 0-25 is shown in figure 1. 


¥ 
0204 


—>x 





‘O1P 





FIGURE 1. Profile of the surface of the standing wave at t= (n+ 4) 7, for n even (solid 
curve) and n odd (broken curve). These curves are based on equation (44) with ¢ = 0-05 
and h = 0:25. 


When the depth / becomes infinite, our results (40) and (42) become 
én = (€+ gis€%) cos x + $e? cos 2x + Ze* cos 3x, (46) 
w = 1—te?, (47) 


These results agree with those of Penney & Price who obtained additional terms 
in this case. 

In calculating the pressure from (9) it is convenient to consider those instants 
when the fluid is at rest. Then in the vertical plane x = 0 and at the time of a 
crest, (9) becomes 


k 63 cosh (y+h) 
= —_ — we 9w~ 8 — 234u~4 81—8w4 . 
pg (P—Po) e* [ 256 (90% "> os] cosh h 
(ORR RA Te _g, cosh 2(y +h) 63 : ee 
_ e| Li? + 3(wp — Wo °) a ae + 556 (1 + 3w@) (27w5 
cosh 3(y + h) 


— 6309 § + 39w9 4 — 5 + 2w§) for -h<y<Mmax- (48) 


cosh 3h 
Upon replacing e by —e in (48), we obtain the pressure in the vertical plane z = 0 
at the time of a trough. Graphs of these pressures are shown in figure 2 for 
€ = 0:05 and h = 0-25. 

When the depth h becomes infinite, our result (48) becomes 


k : 
pg PP) = —y + (e— 336°) ev — fe. (49) 
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This result agrees to the third order with the equation (135) of Penney & Price 
with the exception of the term independent of y. The difference arises from the 
fact that Penney & Price did not take into account the contribution of the zero 
order term in their potential function ®. Had they done this by utilizing their 
equation (47), the correct form of their equation (135) to the fifth order would 
have been 
Aps ora Sas coe 
p = Poly 3A? 4 Att (A—HAP + ERAS) ov 
abe qrAte a0... (; 4, A5) e3u ce (33,45) e>”). 


Our solution describes the standing wave which results from the reflexion of 
a normally incident progressing wave from a wall or breakwater at 2 = 0. It 
also describes the free vibrations of the liquid contained between two vertical 
walls one-half wavelength apart. In either of these cases the pressure given by 
(48) and shown in figure 2 is the pressure on the wall. 
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Ficure 2. Pressure as a function of depth at «= 0 for t= (n+ 4)7, n even (upper curve) 
and n odd (lower curve). These curves are based on equation (48) withe = 0-05 and h = 0-255. 
The horizontal scale is k( p — py)/pg and the vertical scale is y. The dashed parts are sketched 
in to make p—p, have the value zero at the surface y=€p,, and ¥ =€)pin- Equation 
(48) does not yield p—p, = 0 at these points since it is correct only through e°. 


Appendix 
Following Penney & Price (1952), we write 7 and ¢ in a form which obviously 
satisfies (1), (4) and (5). It is 


ie) 
en = > a,(t) cos nz, (1) 
n=1 
io @) 
ep= D> £,(t)e™ cosnz. (2) 
n=—® 
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To satisfy (4) the coefficients f,,(t) and #_,,(t) must obey the relation 
enh B_(t) = e-™ B, (t). (3) 


Upon inserting (1) and (2) into (2) and (3) of § 2, we obtain the following equations 


for the determination of «,, #, and o: 


x ; ? ; l x x 
4=- LY wfi[E(n,s—n)+HE(n,s+n)]—5 & Ss mn 
vt I -— xm x 
xB, 8, |E(m+n,s—m+n)+HE(m+n,s+m—n)] for s=0,1,2,..., (4) 
x x x 
wa,= > nf, [E(n,s—n)+E(n,s+n)]-—5 YX YY mna,,f, 
n x “n om=1 


<|E(n,s—m+n)+ E(n,s+m—n)—E(n,s—m—n)—E(n,s+m+n)] 
for s=0,1,2,.... (8) 


Here a prime denotes differentiation with respect to ¢ and the quantity E(A, s) 
is defined by the equations 


, :. a ; 
E(A, +s) = oF any >N(S); (6) 
ow le 6) aw 
Sy +8) = b pH b» vee Ap Ay by --- ky men—pr eres (7) 
m On op L 
Ap, = Oy (= 25253; <:-): (8) 


Instead of using the iteration method of Penney & Price, we solve (4) and (5) 
by assuming that ~,, £,, and w can be expressed as 


20 


ea, = > a,,€/, (9) 

j=s 
6,= > by, (10) 

j='sl 
o= Yo. (11) 

j=0 

Then H(A, s) can be written as 

E\(A, S$) an Dy E,(A, s)eé, (12) 

j=\|s 


We now insert equations (9) to (12) into (4) and (5) and equate the coefficients 
of the powers of ¢. In this way we obtain, for each pair of integers v > s > 0, the 


equations 
vr v j 
2,=- DY YY ZY of -vlH,-;(n,8—n) + £,_;(n,8+n)] 
n vj=|n| k=0 
le oe 
= LY SZ DZD mnayj by lL, (m+n, s—m+n) 
an vm=0 j=mi=|n| 


+E, ;_(m+n,s+m—n)], (13) 


tv—l 








nts 
the 


13) 
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v-s v v 
Db: WO, A Sy 7 b> = np, |£. j(n,8s—n)+E#,_(n,8+n)] 
0 n=—vj=|n 
] v v v v 
. . . . _ Q ’ - , ip " ; 
9 >> >> p> p> mn mj Pnil Ly i—j\2, 8— M + n)+£, i j(%, 8 +m—*) 
aa vm=VjJ—m 2 n 


—E,,_;_(n,s—m—n)—E,_;_(n,s+m+n)]. (14) 


These equations can be solved successively, starting with s = v = 0. We have 
solved them for all s and v satisfying 3 > v > s > 0. The results are the same 
as those of §3, but the labour is greater. 


The research reported in this paper has been sponsored by the Office of Naval 
Research under Contract No. (285) 45. 
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Exact one-dimensional solutions of the magnetohydrodynamic equations of 
an incompressible fluid are considered. It is shown that one class of plane wave 
solutions of the linearized equations is also a possible class of solutions of the 
general equations including the effect of displacement current. A similar result 
is also established for the solutions for a horizontally stratified fluid. For the 
particular case when the viscosity is equal to the magnetic diffusivity an exact 
solution is obtained for the magnetohydrodynamic Rayleigh problem for a semi- 
infinite plate. It is shown that this solution may be employed directly to give the 
solution for liquids of small, but not necessarily equal, viscosity and magnetic 
diffusivity. 


1, Introduction 

The solution of the linearized equations of magnetohydrodynamics is a topic 
which has received a considerable amount of attention in recent years. The most 
comprehensive and systematic work on this subject is that of Banos (1955a, b), 
who has considered the detailed form of plane wave solutions of the linearized 
equations. It has also been observed by certain authors (e.g. Lundquist 1952) 
that, for an incompressible fluid, the general non-linear equations also possess 
plane wave solutions similar to those of the linearized equations. It thus seems 
of interest to consider other conditions under which the general and linear 
equations of an incompressible fluid possess similar exact solutions. 

In the present paper we examine whether there exist similar solutions of the 
linear and exact equations of motion of an incompressible fluid which are func- 
tions of one Cartesian variable z and time. 

In §2 we consider the particular class of solutions with no velocity component 
in the z-direction. If the displacement current is neglected the exact equations 
are identical with the linearized ones and thus exact solutions of the linear equa- 
tions are also solutions of the complete equations. In particular the forms of the 
plane wave solutions of the linear and exact equations are identical. This fact 
has been observed by Lundquist and others. It is also shown that, without neglec- 
ting the displacement current, there exist solutions of the linear equations which 
are also solutions of the exact equations. For plane waves these are the pressure 
modes defined by Banos (19555) and are characterized by the fact that the as- 
sociated Poynting vector is in the direction of propagation. 

In §3 the particular case when the external magnetic field is in the z direction 
is considered. It is shown that, for a horizontally stratified incompressible fluid, 
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the small amplitude Alfvén wave type solutions obtained by Ferraro (1954) 
represent possible exact solutions of the equations of motion with the effect of 
displacement current included. It is also shown that the equations of motion for 
one-dimensional disturbances in an unbounded fluid are essentially linear and 
capable of exact solution. 

The type of solutions considered in §2 occur in some initial value problems in 
magnetohydrodynamics; one such problem is the Rayleigh problem for an in- 
finite plate when there is a transverse field present perpendicular to the plate. 
The solution of this problem has been considered in detail by Ludford (1959) and 
Chang & Yen (1959). Another class of boundary-value problem satisfying our 
conditions is the problem of the motion of a viscous fluid bounded by an infinite 
oscillating plane, and the solution for an insulating plane has been given by Kaku- 
tani (1959). The magnetohydrodynamic Rayleigh problem for a half-plane is 
another case when the general equations are exactly reducible to a linear form. 
This problem is considered in §4 where it is shown that, when the viscosity is 
equal to the magnetic diffusivity (conductivity x permeability)—!, the boundary- 
value problem reduces to a classical one. It is also shown that the first-order 
solution for small viscosity and diffusivity may be deduced immediately from 
this special case. 


2. General equations 
The motion of an incompressible conducting fluid under the action of a con- 
stant external magnetic field H, is governed by the equations 


ke = = —curlE, (1) 

ct 

] 
(J —p,.V) = E+yv(H,+H), ”) 
dv ° ‘ 
P= = —gradp+p,E+ "J x (Hy+H)+AV°v, (3) 
J = ntti. (4) 
ot 

p, = edivE, (5) 
div v = 0, (6) 
“+ V-grad p = 0, (7) 


In the above equations E and H denote the induced electric and magnetic fields, 
J denotes the current density, v the fluid velocity and p,p and p, the pressure, 
fluid density and charge density, respectively. It will also be assumed that 
(permeability), €¢ (dielectric constant), o (conductivity) and A (coefficient of 
viscosity) are constant. We also have from Maxwell’s equation that 


div H = 0. (8) 
We now investigate the possibility of obtaining exact solutions of the above 


equations in which the field components are functions of one Cartesian variable z 
and of time. Equations (1) and (8) then show that H, is constant; it may thus be 
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absorbed into the z-component of H, and hence there will be no loss of generality 
in taking H, = 0. Equations (6) and (7) have, as one solution, 


v, = 0p ot = dp/oz = 0 


and we consider solutions of equations (1)—(5) satisfying these conditions. Equa- 
tion (7) shows that p will be constant if constant at any particular instant. 
For the first part of the investigation, it will be assumed that the displacement 
current may be neglected (i.e. p, = 9). 
The above assumptions enable (3) to be rewritten 
A A “3. 
OV C : 0*v oH 
fa, > > k (Pp ‘i d/(Hy zi H)?] +i, + pH, ae (9) 
et Oz oz" Cz 
where k is the unit vector in the z-direction. Equations (1), (2), (4) and (8) now 
give 


| | c?H oH ov 
curl J = oa = PG + fH ge = - (10) 
Since », = H, = 0 we have, from (9), that 
5, (P+ (Hy + Hy} = 0. (11) 


H may now be eliminated from (9) and (10) to give 


A 72 27 
(ue -- v= 28 “ead (12) 
ot aodz* C2" 

Equation (12) has been obtained by Ludford (1959) for the particular case when 
v and H, are perpendicular, it has also been derived by Dungey (1958) by linear- 
izing the equations of motion. 

In view of the fact that our assumptions concerning the form of the solutions 
have reduced the equations to a linear form it seems of interest to examine 
whether there exist solutions of the linearized equations which satisfy our as- 
sumptions and which may therefore represent solutions of the general equations. 
Since z is arbitrary it is seen that plane waves satisfy our assumptions, and the 
propagation of small amplitude plane waves has been considered in detail by 
Banos (1955a, b). Banos’s work on magnetohydrodynamic waves is confined to 
the case of an inviscid fluid; ina later investigation (Banos 1956) he has considered 
the propagation of magneto-elastic waves. With suitable changes of notation 
the latter work may be used to obtain the results for the propagation of small 
amplitude plane waves in a viscous fluid. 

Banos has shown that the solution of the linearized equations can be split up 


into two independent classes of solutions which he calls v- and p-modes. The 
v-modes are solutions which are obtained by requiring that the velocity be per- 
pendicular to both k and H,; the pressure associated with such modes is then 
shown to be zero. The p-mode solutions arise when the fluid velocity is in the 
plane of k and H, and perpendicular to k; in general the pressure associated with 
this mode is non-zero. Another difference between the v- and p-modes is that the 
Poynting vector for the v-modes is in the direction of H, whilst that for the 
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p-modes is parallel to k. If the displacement current is neglected the wave- 
numbers of the v- and p-modes are equal. 

We now examine the extent to which Bafios’s results are applicable to plane 
wave solutions of the general equations. If it is assumed that v and H are both 
proportional to exp i(wt + kz), then equation (12) gives 


(iwp + Ak*) (ion + ko) = — PAAR R. (13) 


For A = 0 equation (13) reduces to one obtained by Banos (19555) and for non- 
zero A may be identified with a similar equation derived by Banos (1956) for 
magneto-elastic waves. Thus, if the displacement current is neglected, the wave- 
numbers of the plane wave solutions of the exact and linearized equations are 
equal. 

The existence of v- and p-mode types of solution of the general equations will 
now be considered. It is easily seen that equation (9) may be re-written 


p= = ~te (p+}uH?) +A + nd x Hy. (14) 
é 0z oz? 
Clearly a possible class of solutions of (14) is obtained by setting p = —4uH?; 


v will then be perpendicular to Hy. Hence, since v, = 0, the velocity is per- 
pendicular to both k and Hp. This is precisely the linear v-mode and hence we 
conclude that the exact equations, neglecting displacement currents, possess 
plane wave solutions of the v-mode type. In the exact solution, however, the 
pressure associated with this mode is — 4H? and not zero; in the linear theory 
terms of order H? are neglected and this would give p = 0. Equation (13) shows 
that there will be two possible v-mode solutions, one of which is essentially a 
damped Alfvén wave and the other a highly attenuated one which vanishes for 
zero viscosity or infinite conductivity. v will be of the form Ak x H, exp (wt + kz) 
and H and E may be obtained from equations (1) and (10). 

It would be logical at this point to consider the existence of p-mode type solu- 
tions. It will, however, be shown that in this case it is not necessary to neglect 
the displacement current and we therefore consider conditions under which the 
effect of the displacement current may be treated simply. The non-linear effect 
of the displacement current is due to the charge density term and thus, in 
order to seek simple solutions, we shall consider only those with div E =0. 
Since our main interest lies in solutions of the plane wave type it may be assumed 
that EF, = 0. If the displacement current is included then (10) becomes 


1 o? c?H oH : (15) 


C 
wi —~ 6 = — A, = 
Co oz C te ot “@ 


~ 


Equation (15) shows that v and H are parallel for plane waves, and equations 
(2) and (3) thus become 


l 
= E+ 6 
= E+yvxH, (16) 
OV “ ov oE i 
px = —grad(p+4uH?)+p"J x Hy)+A~, — fe x H. (17) 
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The only non-vanishing component of the last term on the right-hand side of 


(17) is that in the z-direction and hence this term may be written as —¢ grad yy, 
where vy is a scalar depending on z and ¢t, and hence 


o = —grad(p+suH?+eyr)+A tn x Hh). (18) 
Equations (2) and (4) give 
~J.Hy = E.Hy, (19) 
(1- ie ) J.Hy = Hy. curl H. (20) 
From (1) and (19) 
iF (3 } a aiba\S-He = ~ “it 5 (div). (21) 


Thus, since div E = 0, J.H, = E.H, = 0 is a possible solution and E will have 
the form ¢k x H, where ¢ is a scalar. From (16) 


(p = = a) curl v = «curl (J x H)), (22) 


02° 
and eliminating v between (15) and (21) gives 


eJ 0E)\] 2, 2 





C ( jl . 
SS — curl | curl ioe ll - Ay, =, =, curl (J x H,). (23) 
where v is the kinematic viscosity. 
Finally, from (1), (4) and (22), 
0 e\( 1 (2a 02\ 0)0?d HG, co? (0? e2 ; 
(a? =) | on (52-5) Tal de poe (s2- aa)? — (24) 


If the displacement current is not neglected then (9) has the form 


pW = — KS [p+ Ju(Hy +H)? + ey] + ek He — HE) 
+ 6tHp, [ke x (k x H,)] i + Aa (25) 
Vv is perpendicular to k and thus from (24) 
ov e2y rn) cH 
“a = Ax + eH, [K x (k x H,)] + Hy, ae? (26) 
C od 
C [p+ Yye(Hy +H)? + ey] = el He — HEY. (27) 
C2 Co 


The velocity and magnetic field components are now determined in terms of ¢ 
from (25) and (26) and the pressure is then given by (27). 
If it is assumed that ¢ is proportional to exp i(wt + kz) then (24) becomes 


(1 
lou 


| Hy, 


(iw + vk?) (k? — enw?) + tw) = —= (enw? — k?). (28) 
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For v = 0, (28) reduces to the dispersion equation obtained by Bajos for the 
p-modes of an inviscid fluid and with suitable changes of notation is identical 
with the corresponding equation obtained for magneto-elastic waves. For 
plane waves v and H are parallel and thus equation (26) shows that v is parallel 
to k x (kK x H,); this is precisely Baiios’s definition of the p-mode solution and the 
other field components are identical with the ones obtained by Bajos for the 
linearized equations. Thus the p-mode solutions of the linearized equations are 
possible solutions of the general magnetohydrodynamic equations. If H, and 
k are parallel, the appropriate form for E is then ¢a x H, and it may then be 
verified that in this case the p- and v-modes are indistinguishable. For H, and k 
parallel it is possible to obtain some additional exact solutions of the equations 
and this point will be considered in more detail in the following section. 

The fact that the present method of generating p-mode solutions for plane 
waves is identical with Banos’s is a particular example of a general result 
(Williams 1960) that the general solution of the linearized equations can be ex- 
pressed in terms of two independent classes of solutions. One class is generated 
by a stream function and the other by a one component electric vector potential. 
For plane waves the first class becomes the v-mode solution and the other the 
p-mode solution. 


3. Wave fronts perpendicular to external fields 

The problem considered in the first part of the present section is of a more 
general nature than that of the previous sections in that we investigate the 
deviations from the uniform state E = H = 0 of a horizontally stratified liquid. 
[t is also assumed that the permanent magnetic field is vertical. The propagation 
of small amplitude disturbances in an incompressible fluid of this type has been 
examined by Ferraro (1954); the corresponding problem for a compressible fluid 
has been solved bv Ferraro & Plumpton (1958). Ferraro has shown that for an 
incompressible fluid there will be propagated small amplitude disturbances 
which are essentially Alfvén waves with a variable velocity. It will now be shown 
that the exact equations also possess solutions of this type. 

The undisturbed pressure and density will be denoted by p, and pg, respectively, 
and the deviations from these quantities by p, and p,. Equation (3) now becomes 


y 


dv , 
(Po +P) > — grad (py+p,)+p, div E + uJ x (Hy +H) — (p9+ p1) 9k, 


(29) 
where g is the acceleration due to gravity. The static equilibrium condition is 
Apo ‘ 

dz = —Pod; (30) 


In order to simplify the analysis it has been assumed that the fluid is inviscid; 
this assumption is not actually essential and the analysis can in fact be carried 
out for a varying coefficient of viscosity. One solution of equations (6)—(8) is 
v, = H, = p, = 0, and we consider solutions of the equations which satisfy this 
condition. 
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We now attempt to obtain a particular solution of the equations with v and H 
parallel; this assumption is made because the solution obtained by Ferraro 
satisfies it and also it enables the effect of displacement currents to be included. 
Equations (2) and (4) now show that 


H 


€ ¢ . ~ 
(1 + an) iv E = 0, (31) 


and a possible class of solutions of (31) will be div E = 0 = E.. From (1), (2) and 


(4) it may be shown that 





a 1/82 — a\ a gt 

é ey! = - ry 32 
MH ata asa) [J oe et —5| V x Hy, (32) 

and (29) and (32) then give 
0 1/2 o2 (oe C2 
les = ‘a) | = — Hi} (= 3 Ell ) V 
- @ 1/2 02 

- —- — P 3¢ 

= lH shes cu 5a) | grad P (33) 


where P= p+4uH?+ey and y is defined as in (18). Since v, = 0 and p is a 
function of z and t only we have that the right-hand side of (33) vanishes and hence 


2 (0 02 0 \ 
Posie - (= — El eo (fo — +°H3)| v= 0. (34) 
v will have the form da x k, where ¢ is a scalar solution of (34) and a an arbi- 
trary constant vector; the other components may then be expressed in terms of ¢. 

For ¢ = a! = 0 equation (34) reduces to the Alfvén wave equation with a 
variable velocity (“H?/p,) and is the equation obtained by Ferraro by linear- 
izing the equations. Thus we see that for this class of problem, also, there exist 
solutions of the linearized equations which also represent possible solutions of 
the general equations. 

The above analysis has been concerned more with particular forced solutions 
of the general equations than with the solution of particular boundary-value 
problems. It is also of interest, however, to examine the possibility of obtaining 
exact solutions for initial value problems where the wave fronts are perpen- 
dicular to the external magnetic field. The fluid is assumed to occupy an un- 
bounded region of space and the initial density is assumed to be uniform; equation 
(7) then shows that the density will always be uniform. Equation (3) becomes 


C C oH, ev 
P| + v(t) ly — —grad [p+ 4u(H, + H)?) + Ah. = +X 3 ; (35) 


ct - FA a 
From equations (6) and (8) it is seen that the gradient term in (35) vanishes in an 
unbounded region and hence that év,/ct = 0 and v, is thus constant. From (2), 


4) and (8 
ee cH dH 


= — 


+H —U,=—)5 (36) 
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where the displacement current has been neglected in deriving (36). Equation 
(35) now becomes 


(c Cc 0? | oH _ 
{— +y,— —yp-_\v = nH, —. 37 
F \ct” 702 ez3 i 7) 
Eliminating H from (36) and (37) gives 
0 C it @ ¢ o? | > O-V ; 
i= +0, —Vas) (Ue ta a VC = AG = 38) 
P\dt* “*2z O22) | *ez ct aomdz} Me 578 = 


Equation (38) is non-linear but the non-linear effects are due to the constant 
term v, and thus the equation is amenable to exact solution. For v = 0 this 
equation (without the v,) was obtained by Roberts (1955) by linearizing the equa- 
tions of motion and employed by him to solve an initial value problem with 
v, = 0. The above analysis thus shows that Roberts’s solution in fact represents 
an exact solution of the magnetohydrodynamic equations neglecting displace- 
ment current. 

In the following section the above analysis will be applied to obtain the 
magnetohydrodynamic analogues of some exact solutions of the Navier-Stokes 
equations. 


4. Viscous flow problems 


One of the classical exact solutions of the Navier-Stokes equations is that 
for the Rayleigh problem. This problem is essentially the solution of the equations 
of motion for a viscous incompressible fluid outside an infinite flat plate which is 
suddenly moved parallel to its length with a velocity U. The corresponding 
problem when the plate is perfectly conducting and there is a uniform magnetic 
field present perpendicular to the plate has been solved by Ludford (1959) and 
Chang & Yen (1959). The plate is assumed to occupy the plane z = 0 and to be 
moved parallel to itself with velocity U in the x-direction. Clearly this problem 
is of the type considered in §3 and v, will satisfy (12). The solution obtained by 
the above authors, by using Laplace transformations, is rather complicated 
but takes on a considerably simpler form for v = (uo)~!. It is also of interest 
to note that, for a perfectly conducting fluid, the solution may be obtained 
from that of a completely different type of initial value problem solved by 
Roberts (1955). 

It is of mathematical interest to note that equation (38) may be used to solve 
the magnetohydrodynamic Rayleigh problem with constant suction on the plate. 
This problem has been solved in the absence of a magnetic field by Hasimoto 
(1956). 

The Rayleigh problem for a semi-infinite plate may also be reduced to the 
solution of a simple linear equation. We assume that the plate occupies the 
region z = 0,y > 0, for all x; it is then easily verified that a solution of the equa- 


tions of motion is possible with v, = v, = 0, and v,, a function of y, z and f, 
J zZ Ir e 
satisfying 
- 7 72 
C c 070 , 
(- — V8) (< — V2)», = V2=Z. (39) 
ct } \ot \ oz" 
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9 = (uo)~!, V2 is the two-dimensional Laplacian operator and V the Alfvén velo- 
city. The only component of H is H, and equations (9) and (10) become 


(nv?-<) H, + Hy? = 0. (41) 


~ 
~ 


The electromagnetic boundary conditions are that the tangential components 
of E and H and the normal components of “H are continuous at the plate. It 
can be shown that the boundary conditions on the plate reduce to v, = U and 
n OH,/Cz = 0. For infinitely conducting fluid the second boundary condition is 
satisfied identically and in this case the only condition is v, = U. The general 
boundary-value problem is rather complicated; the Laplace transform method 
is clearly the appropriate one to employ. The problem is then clearly of the 
Wiener—Hopf type, the transform variable p will, however, occur in the Wiener- 
Hopf factorization and the resulting general solution will be extremely com- 
plicated. It will, however, be shown that for v = 7 = k,, the boundary-value 
problem reduces to one of a type already solved. 

In this case equation (39) becomes 

(4,V?-—-Vo) (2-540; 
\ ct U C 


2 
Thus, from (42), we see that 

v, = eH f(y, z,t) + e-"2/*h gly, z, t) (43) 
where f and g are independent solutions of 


0 V2 
v — ee ain = ), 
(! Ct s sige 


Equation (40) now gives 
Ay 5 yaar risa 
H, = — y {fer 22h —geV 22k}, (44) 
Equations (43) and (44) show that the boundary conditions on the half-plane 
become ‘ y 
» y CO y e 
f+g=U and ~(f—g)t+5 (+9) = 9. 
If 5 and 1 are defined by 
@g=ft+g and w=q-f, 
then the boundary conditions on ¢ and y on the half plane are 


ee VWlao — WV ik 
@?=U and dy/oz = UV/2k,. 


The boundary-value problems for ¢ and y are of a well-known type and may be 
deduced from the solution of the ordinary Rayleigh problem for a half-plane 
solved by Howarth (1950). The best method of effecting a solution is by means of 
the Laplace transform. If w(y,z,p) denotes the Laplace transform of u(y, z, f), 


then ™ 
“= [ e~Pty dt. 
~0 
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We thus require solutions ¢, y of 


9 : ie 7 D 
{8 (r+32)) u= V0, (45) 


such that 6 = U/p and éy/¢z = UV /2k,p on the half-plane. For V = 0 Howarth 
has obtained a solution ¢,, of (45) such that d,, = U/p on the half plane. This 
solution is defined by 


U 
$n = 55 {2cosh qz—e-*erf q}(y —£) — ew erf q4(y + 8)}, (46) 


where q? = p/k, and y = £?—7? and z = 2£y. 
A solution to the present boundary-value problem for ¢ is thus 


V2\ 1 . 
5=(p ky +) ouly, z,pt+ V 2/4k,) 
and hence 4 


“t 
p=er™ idnlye d+ se | eF?witk, 5 (y, z, w) dw. (47) 
1/6 


By a method similar to that of Howarth’s it is seen that 


y= iar {2 sinh fz —e* erf £3(y +) +e" erf £4(y—£)}, (48) 
where k, f? = p+ V?/4k,. 


The inverses of d and y may be expressed in terms of integrals of functions 
obtained by Howarth, which are themselves infinite integrals involving modified 
Hankel functions. The actual forms of ¢ and y are extremely complicated and 
will therefore not be considered. It is of some interest to examine the behaviour 
of the solution for large and small values of t; for these cases it is possible to 
obtain reasonable simple forms for the solution. A quantity which is of physical 
interest is the skin friction 7 at the wall and we shall consider the limiting forms 
of 7, keeping y fixed, for large and small t. We have that 

Ans a Ff 
” r= {uth He+pok |) = = pak E+ SE “| 


Cz lwall 


On the wall 


Ad ; (Od a d 
Hi = e-V/Ak, ( =) + BS L(y, 2, V /2k | “ae an e-Vtw/4k, ( Ft) dw 
Ce Of } wall | 6 wall Cz / wall 
(49) 
ad : 
For large t, =” ~ 0-4604U (k, ty?)-* 
7 0 wall 


and thus 


OPA 99 var, | V aes — i) r 
al <p ~ 09208 (k, ty2)-t e “tae erf ma y (50) 


wall 
One clearly expects that 0¢/0z will tend to cai as V + 0 but equation (50) 
shows that the asymptotic form of 0¢/éz is twice that of 0¢,,/0z in the limiting 
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case of small V. This apparent discrepancy is due to the fact that for V small 
the integral in (49) may not be estimated by normal methods and a method 
developed by Clemmow (1950) must be employed. It can be shown that the 





appropriate form of the first term in (50) is 


1 
0-4604 - Vt w 
| eV ak, +( ) [ © de, (51) _ 
(i, ty”)! | 4h, VU 4k, we } ¢ 
Equation (51) exhibits the appropriate limiting behaviour for V = 0 and the 
integral is an incomplete gamma function. For small ¢ it can be shown that ‘ 
04 Vt : 
eg xr ; | 
/ wall 5 eae 4h, 
€ 
It is easily shown by transform techniques that C 
y= -Vert [oes O(t-t) 
V 
UVt 
for large ¢, and y= + O(t) ( 
1k, 
for small ¢t. Thus we finally obtain 
| 2! f Vy / it Vy\ 0-9208 li ‘ 
T= —Pok er air if exp { - = — pexp ( - ; as > 
Pom” hy 2k, | A 2ky E 2k,) (y®k, t)t J 4k,}} | | 
Ink, U V% 
er) ae ESM, | 
Jmkt | 4k) 
(52) 
For small values of ¢t the skin friction is that obtained by Howarth. The first two T 
terms in (52) represent a residual skin friction; there are two contributions, one J 
from the residual magnetic stress and the other from a residual viscous stress. | 
A similar state of affairs occurs in Ludford’s work, the only difference being the | a) 
absence of a residual viscous stress. The third term in (52) is 2exp (— V*t/4k,) 
times the corresponding skin friction obtained by Howarth, and for small V this z 
term should be replaced by equation (51). The difficulties encountered in the b 
expansion of 7 for large values of t may be clarified by writing 7 as py) U(k,/t)? F(R, S) d 
where R? = y?/k,t, S = V*t/k,. It is now seen that (52) gives the form of 7 for SI 


small R assuming S large, whilst equation (51) gives the modification necessary | 
for finite values of S, still assuming FR large. A similar situation occurs in Lud- 
ford’s work and since there appears to be a slight discrepancy in his work we shall fc 
consider this point very briefly. Ludford obtains 


aT ne ae 
= Verf fe exp (— i. a. 
pl Vatu" a at (Jn +»)? 


and thus for large ¢ 
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The skin friction is thus 

2 + Ji exp — Vt/(./y +./v)? 
times the corresponding result obtained by Rayleigh and again a discrepancy 
occurs in the limiting case of small V. This discrepancy is due to the fact that for 
small V it is not possible to expand the error function asymptotically. Ludford 
has obtained — 1 instead of {2 + ./7/v}. 

Clearly for 7 and vy small we can obtain a solution by perturbation methods 
expanding in powers of 7 — v. The process will be fairly complicated but it will be 
shown that a solution neglecting squares of (7 —v) may be obtained immediately 
from the solution for v = 7. If k, is now defined by $(7 +) then (39) shows that 
equation (42) will still hold if terms of order (7 — v)? are neglected. With this new 
definition of k, we define V’. by 


Vio = eV etki ft e Veh g, 


where f and g satisfy the same boundary conditions as above. We also have from 
(4) that as - 
y2 _ +H, (i v2 4 V+ 4H(v—9) V2V2 = 0. (53) 
5 


Thus, neglecting terms of order (7 — v)? 


,0H HV m 
ye (+914 Hox (fa) + AV Hy y=) 
0 V r 
«{gU-9)+5¢ +9) +4 H, oY —1)) =, (f+g) = 


The conditions on f and g show that, to order (v — 9)?, 0H//éz = 0 on the wall. Thus 
V; represents a solution of the problem neglecting terms of order (7 — v)?. 

For 7 = 0 the solution may be obtained from Howarth’s by a suitable variation 
of the parameters. 

Another viscous flow problem capable of exact solutions is that when the plane 

= 0 oscillates parallel to itself. For an insulating wall a solution has been given 
by Kakutani (1959). The solution for an oscillating half-plane may also be de- 
duced immediately from the Laplace transform of the above solution for impul- 
sive motion. 


The author is indebted to a referee for some detailed comments on an earlier 
form of this paper and for some additional references. 
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Starting with expressions for viscous stress, heat flux and diffusion flux, we 
formulate a continuum theory for steady flow of a binary mixture of chemically 
inert perfect gases through a normal shock of arbitrary strength. For shocks of 
vanishing strength, a solution by series expansion in Grad’s (1952) shock-strength 
parameter gives a result essentially the same as found previously by Dyakov 
(1954). For stronger shocks a straightforward numerical integration, quite 
analogous to that useful in the simpler pure-gas problem, is laid out. 

The resulting problem has eight parameters: shock strength, ratio of specific 
heats, ratio of bulk viscosity to shear viscosity, Prandtl number, Schmidt number, 
thermal diffusion factor, molecular mass ratio, and initial mixture concentration. 
A dozen examples were worked out on a simple desk calculator to exhibit the 
influence of some of these parameters. They involve the gas pairs argon*— 
argon*®, argon—neon, argon—helium, and xenon-helium. 

In discussing the results, special attention is paid to the degree of success with 
which the weak-shock theory may be extrapolated to arbitrary shock strength, 
and to the question of the accuracy of the Navier-Stokes approximation for a 
mixture of gases of very different molecular weights. 


1. Introduction 

A shock wave in a pure gas is broadened into a continuous transition zone by 
viscosity and thermal conduction, in the view of the continuum theory of gases. 
From this same viewpoint, a shock in a gas mixture is further broadened by 
dissipative species-diffusion processes, and a partial separation of the mixture is 
caused by the gradients within the shock. 

Previous theoretical studies of shock waves and sound waves in binary gas 
mixtures have been made, notably by Cowling (1942), Dyakov (1954), and 
Kohler (1949), and the essential results for very weak waves have found their 
way into important reference works by Herzfeld (1955) and Hirschfelder, Curtiss 
& Bird (1954). Interesting and somewhat related papers on shocks in dusty or 
ionized gases have recently been contributed by Carrier (1958), Jukes (1957), and 
Tidman (1958). The present paper is primarily an extension and elaboration of 
the works of Cowling and Dyakov. Whereas Cowling neglected viscosity, thermal 
conduction and thermal diffusion in order to gain a preliminary insight into 
diffusion effects in shocks, and Dyakov limited his treatment to waves of sufficient 
weakness so that the separation of the mixture could be neglected as a higher order 
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effect at the first stage of calculation, we make a direct frontal attack by numerical 
integration, to find the structure of a shock of arbitrary strength in a binary 
mixture of chemically inert perfect gases, including effects of viscosity, thermal 
conduction, ordinary-, baro-, and thermal-diffusion. 

Of course, even accurate numerical integrations do not make the present theory 
(based on the Navier-Stokes approximations to viscous stress, heat flux, and 
diffusion flux) valid for arbitrary shock strength. Indeed one might intuitively 
expect that the maximum shock strength for which these calculations will give 
realistic predictions would be reduced as the molecular mass ratio of the mixture 
increases, due to the difficulty of maintaining thermal equilibrium between the 
light and the heavy constituents. Thus experiment may prove someday that the 
more ‘interesting’ features of the present results, which mostly appear when the 
mass ratio is fairly large and the shock strength is finite, are not realistic at all. 
Nevertheless, it has been interesting and instructive to the author to see just what 
the Navier-Stokes approximation says about this rather complicated problem. 


2. Diffusion velocity and conservation of mass 

We take as our starting point the results of Chapman and Enskog for the 
diffusion velocities of a binary mixture of perfect gases, obtained by kinetic 
theory. Equations 8.4, 7 and 8.3, 7 of Chapman & Cowling (1939) may be com- 
bined with the definitions of the total density and mass velocity of the mixture, 


P= Pi +P» (1) 
PU = PyUyt+PoUs, (2) 
to give, in our present nomenclature, 


_ pM, My, {, M,-M, d dm] a) ; 
1 ee ee ee ee al 


Subscript 1 refers to the heavier molecular species and / is the mole-fraction of that 


f =p, M/eM,. (4) 


species 


The mean molecular mass M of the mixture depends on f, being given by 
M = fM,+(1—f) dh. (5) 


The binary diffusion coefficient D,, and the thermal-diffusion factor « are positive 
quantities, except in some exceptional gas pairs for which « may be slightly 
negative. 

The qualitative nature of the separation of species within a normal shock can 
be seen directly from (3). If we view the shock as fixed in space, with steady flow 
from left to right so that uw > 0, we then expect dp/dx > 0 and dT /dx > 0 within 
the shock (unless the situation varies drastically from that in a pure gas). The 
driving forces tending to separate the mixture are d(Inp)/dx and d(InT)/dz, 
leading respectively to baro-diffusion and thermal diffusion. From (2) we see 


that baro-diffusion speeds up the heavier component, M,—M, > 0, and slows down 
the lighter component relative to the mass velocity of the mixture. Ife > 0, which 
we would normally expect for J/, not too close to M, (ef. Chapman & Cowling 
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1952, p. 254), then thermal diffusion slows down the heavier component and speeds 
up the lighter. (This is consistent with the sign convention for «, namely, that it is 
positive when the heavier molecules seek the colder regions.) Thus under most 
cirumstances thermal diffusion will partially counteract baro-diffusion within 
a shock wave. This same conclusion was reached by Dyakov (1954) and, in the 
related acoustic problem, by Kohler (1949). Of course, once any separation has 
taken place, ordinary diffusion, with driving potential df/dz, will tend to homo- 
genize the mixture. 

In the absence of chemical reactions each component of the mixture is subject 
to a simple continuity equation, 


PU, = Mm, (constant), (6) 
PoUy = My (constant). (7) 

Equation (2) gives 
pu = m,+m, = m (constant) (8) 


as the continuity equation of the mixture. These equations make it clear that 
slowing down a component relative to the mixture velocity increases the con- 
centration of that component relative to its value outside the shock wave. Thus, 
if baro-diffusion dominates thermal diffusion, the lighter molecules will be con- 
centrated in the shock. This is usually the case. The degree of concentration would 
be expected to be greatest about where d(In p)/dz is maximum. For very weak 
shocks this proves to be at the centre of the shock, and the concentration profile is 
symmetric about this point (Dyakov 1954). For stronger shocks, the maximum 
d(In p)/dx will be upstream of the maximum of dp/dz, so we expect a peak con- 
centration of the light species upstream of the centre of the shock. This upstream 
shift of the separation maximum is reinforced by ordinary diffusion, which 
accelerates the heavier molecules relative to the mixture when they are upstream 
of the separation maximum and decelerates them after they pass it. 

Note that equations (3), (6) and (8) do not suggest that the relative scarcity of 
heavy molecules in the upstream portion of the shock is compensated by an equal 
enrichment of this species at some other position in —0« < 2 < 0, unless dp/dx 
should ever change sign. It does not in any of the examples worked out here. The 
tempting intuitive notion that conservation of mass prohibits a net deficiency of 
either species relative to the proportions of the mixture at infinity is not applicable 
to steady-state analysis, and usually arises when we include in our thinking the 
essentially transient idea of shock generation in a mixture of initially homogeneous 
composition. 

An accurate analysis of transient shock formation in a mixture is outside our 
present scope. However, a qualitative description is suggested now in hopes of 
clarifying the relation between the following strictly steady-state results and the 
corresponding ‘long-time’ behaviour of a particular shock-formation problem. 
Suppose that a piston is impulsively accelerated to constant velocity in a tube 
containing a mixture initially homogeneous and at rest. The author supposes that 
in the initial instant (less than a mean molecular collision time) of that piston 
motion, the light molecules struck by the piston fly away from it more rapidly 
than their heavy neighbours. In the course of subsequent collisions, which 
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transmit the impulse of the piston out into the gas, this behaviour is reproduced 
at the front of the column of gas which has started to move with a mean velocity 
equal to that of the piston, and becomes a characteristic feature of the structure 
of this front (the shock wave). After a sufficient number of collisions, the motion 
and structure of this shock front become locally stationary if viewed by an 
observer moving with the shock speed. 

In the meantime, the heavy molecules which move away slowly from the 
piston during the initial instant produce a local surplus of this species in the im- 
mediate vicinity of the piston. This surplus is embedded in a region of gas which 
soon has the mean motion of the piston and redistributes itself by ordinary 
diffusion. The region throughout which this original surplus is subsequently 
spread thus grows out from the piston face in proportion to , (time), while the 
shock front moves away from the piston in proportion to time. Between the 
shock and the region of redistribution of the heavy molecules which are ‘left 
behind’ during shock formation, there appears an ever-lengthening column of 
homogeneous mixture having the velocity of the piston and the composition of 
the undisturbed gas. In the steady-flow analysis of the shock structure, the 
properties of this column are those supposed to exist at z = +00, so that the 
residual effects of shock formation, including the local enrichment of heavy 
molecules which compensates for the local depletion in the shock, are specifically 
excluded from the steady-state theory. 


3. Conservation of momentum 


The diffusion equation (1) results from the momentum conservation principle 
applied to the relative motion of the constituents of the mixture. The momentum 
equation for the mixture moving as a whole is 


du 
dx 


pu* + p— (3H +k) = P = py) + Po- (9) 
which differs from that for a pure gas only through the dependence of # and « on 
fin the present case. In fact, virtually nothing is known about the bulk viscosity 
coefficient « for gas mixtures, since this additional dissipative mechanism is hard 
to separate experimentally from ordinary viscosity, heat conduction, and dif- 
fusion. We carry the symbol « in this analysis just to be open-minded, although 
the subsequent calculation will assume x/ to be constant, with no particular 
physical justification. 


4. Conservation of energy 
The differential equation of conservation of energy for the mixture may be 
integrated to give 


l 
: +q, = 4Q = m(hy+ 4u?). (10) 


9 
pu(h+ uy?) = u(Sye 1 x)=. 


The specific enthalpy in (10) is given by the rule for mixtures of chemically 
inert perfect gases 


ph = pyhy + pahy = (PyCp1 + Popo) T. (11) 
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In (11) explicit use is made of the important assumption that both components of 
the mixture share the same temperature. 

The heat flux q, in (10) contains two terms in addition to the ordinary Fourier 
conduction, one a simple convection of partial enthalpies at the diffusion velocity, 
the other due to the diffusion thermo-effect (Chapman & Cowling 1952, equa- 
tion 8.41, 3), 

dT M? p 
= —-A— +p,(u,—u)|h, -A. +=, a). 12 
Ve dx Py( 1 ( 1 2 M, M, p ( ) 

When (12) is substituted in (10), it proves convenient to group together all the 

enthalpy terms. One easily shows that 


puh + py(u, — u) (hy —hg) = mh, + mah, 


— , ‘yy 
= (MyCpy + MgCyg) 7 


me, T, 


where c,, denotes not the local and variable specific heat of the mixture, but the 
specific heat of the mixture outside the shock (hence a constant for any given 
shock problem). The last step above follows directly from the basic definition of 
specific heat of an inert perfect gas mixture 


Pp a Pi Co + P2ly2 


and the observation that outside the shock u, = u, = u. 
Equation (10) is thus rewritten as 


du _,dT | p,M* 
dx dx pM,M, 


‘yy 


pu(su +c, 7) —u(gut+k) (u,—u) pa = 3Q. (13) 
5. Equation of state 
If we define a gas constant for the mixture outside the shock wave by 


— { ‘jy 
R = polPoT, 
then inside the shock, where the mean molecular mass is variable, we have 


p= a PRT. (14) 


6. Dimensionless variables and boundary conditions 

Equations (1) through (7), (9), (13) and (14) give ten independent relations 
for the ten unknowns /, pj, Po, U, U4, U2, f, M, p, and T’. They can be reduced to 
three coupled first-order differential equations for wu, 7’ and f. These, together with 
a sufficiency of transport property data and the boundary conditions 


du aT df 


= = =Q at r= +0 15 
dx dx dz = se aa 


define the problem. 
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We introduce dimensionless variables and parameters as follows: 


v = mu/P (velocity), 

tT = m*RT'/P? (temperature), 

a = mQ/P? (shock strength), 

y = ¢,/(c, —) (specific heats ratio), 

4 kK 

3 Ps 

§= M,/M, (molecular mass ratio), 

f,= mole fraction of heavy molecules outside the shock, 


and the reference length L = p/m. 


The three equations for v, 7 and f can be arranged as autonomous equations for 
dv/dx, dr/dx and df/dx as follows 


FAs _ l Le (9-1) fy7 a). 
1+(0-l)fv | 


dx X| 
dr y-—1(uc, = - Y 1+(A-1)fy - } a 
Yde = By zr) ee" ~at an ean ey OO 


(16) 


1 e- (fo—f) [1+ (O- DFP f—f) 1+ (O-Df] 
“de —\pDyp) (1+ (P-L +(0-Vfol - 1+(6?—1)f 
x (0-1) © (in v) —[(0—1) —af{1 + (6 gas (n7)}. (18) 


7. Method of solution 


The method used to solve (16)—(18) numerically is quite commonplace in the 
shock-wave literature (e.g. Grad 1952, Gilbarg & Paolucci 1953, or Talbot & 
Sherman 1959) and will only be described in principle here.* 

The independent variable x can be easily eliminated by dividing (17) and (18) 
by (16). This yields two ‘direction field’ equations for d7/dv and df/dv, for which 
we seek a particular solution joining the upstream and downstream states of the 
gas. The co-ordinates of these states (v,7, fata = +00) are easily found by setting 
dv/dx = dr/dx = df/dx = 0 in (16), (17) and (18). 

Since these state points are obviously singular points of the direction-field 
equations, a preliminary investigation of the direction field in their vicinity must 
be made to determine the end-slopes of the desired solution curve and the correct 
point from which to start the numerical integration. This analysis shows in the 
present case that the downstream singular point is almost invariably a saddle 
point and the upstream singular point a node, just as they are for the Navier-— 
Stokes equations in a pure gas. The integration must be started from the down- 
stream singular point and this guarantees the existence of a unique solution. 

After the solution curve is found in v-r-f space, the spatial distribution of v is 
found by a numerical quadrature of (16). 


* Readers interested in these details are invited to communicate directly with the author. 
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8. Concentration dependence of transport properties 


The numerical integration of the direction-field equations and the quadrature 
of (16) require a knowledge of the concentration and temperature dependence of 
viscosity , thermal diffusion factor «, and the ratios ~/A and “/pD,,. Numerical 
data could of course be used directly, but for the present sample calculations 
simple analytic approximations were introduced after a study of data for argon— 
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FIGURE 1. Properties of argon—helium mixtures. 


helium mixtures (cf. figure 1; Chapman & Cowling 1952, pp. 232, 243, 248; 
Hirschfelder et al. 1954, p. 584). These approximations provide a local fit to the 
empirical transport property behaviour in the vicinity of the known mixture 
state ahead of the shock. They are 


ye oc T3, independent of f, 

pA oc M, independent of 7’, 
1/pD,. «< 1/M, independent of 7, 
a o 1/M, independent of 7’. 


While some of these are certainly inaccurate if applied to large variations of f 
and 7’, they are entirely adequate for the purpose at hand, since f actually does 
not vary widely in a given shock, and the assumed temperature dependence is 
fairly realistic except in the case of x. Thermal diffusion is fortunately a sufficiently 
small effect so that sizeable errors in « produce much reduced errors in the shock 
profiles. 
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If we introduce the Prandtl number, Schmidt number, and thermal-diffusion 
factor of the mixture ahead of the shock as basic parameters of the shock-structure 
problem, the approximations above state that 


pc, |A = Pr (M/M) = Pr {1+ (0-1) f)/(1+ (8-1) fol 
(remember that c,, already refers to the mixture outside the shock), 
h/pDy». = Sc(M,/M) 
and a = a,(M,/M). 


9. Series solution for weak shocks 

In the limit of vanishing shock strength, .@ —> 1, a series-expansion solution of 
(16), (17) and (18) is appropriate, and the first terms have been found by Dyakov 
(1954), albeit in a nomenclature quite different from ours. We proceed a little 
differently here, expanding in terms of the shock-strength parameter proposed by 
Grad (1952) instead of the pressure ratio, Dyakov’s choice. One hopes that the 
former expansion parameter, which remains finite as .M — 00, would lead to 
quicker convergence, or at least better approximation by the first term of the 
series, at moderate Mach numbers. In fact, this seems to be the case. 

Grad’s parameter € is introduced by 


1 /y-—1\( 2y . 19 
ti — él, 
: (ye) (y-a+ i 


where the + sign gives vo(z = —o) and the — sign v, (2 = +00). It is related to 
our earlier parameter a by 


; l/y-1\[/ 2vy \?_ , 20 
a= —€* 2 
iF) Fost ey} it 


and to the shock Mach number by 


Dae fee 
jae a4 (21) 
y—-1l\y 4*+1 


Next we define new dimensionless velocity and temperature variables, which 
vary from — 1 to + 1 through the shock, whatever the value of €. These are w and tf, 


/ 


given by 


_ &fy~i 2y : is 
= 305)" @ 

l — 4y . | 
= : 2et —e?}. 2¢ 
; roe \y—-1e FY ne 


A new concentration variable may be introduced by 
24 ¢ 
fo—f = €¢. (24) 
We now expand w, t and ¢ in €¢, i.e., 


w= w+ ew) + ew + ..., 


t= £0) 4 ef) 4 | 2¢(2) 4 eee (25) 


”) —_ Pp + EG + FG + dd 
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and substitute these series into (16), (17) and (18). Grouping terms of like order 
in €, and solving the resulting simple differential equation for w, we get 





y?-1 pur) : 
x = — tanh | > i. i 
u tanh | ; Yat, Fl TSA) W,-M, y-1_\)| x i (26) 
| ‘ Pr Se (= M, yo %o 
{ — — we), (27) 


go = Yt! \(y—1) ” fol (1—fo) 





8y Sc 
| — M, a ‘- | 
' , -\ f(1—(w)2 28 
[fete Yfo(1 —fo) (44 —- +L, en. )| (1 — (w®)?} (28) 
a Pr Se M y *0 


for the lowest-order non-vanishing terms in w, t and ¢. The corresponding result 
for shock thickness based on maximum slope, 


bn, = 2/|dw/da| max> (29) 
pid, __ _8y _ (1\(y,%—1, Woll—fo) (Mh—-M,_y—1, \" is 
He “sale Pr Se | My Y 40) | ai 


When our nomenclature is brought into agreement with that of Dyakov (1954), 
these last four results agree perfectly with his. The quantity in braces in (30) 
is also identical with the corresponding factor in the sound-attenuation coefficient, 
found by Kohler (1949), and cited by Hirschfelder et al. (1954, p. 732, equa- 
tions 11.4-20). It may be noted in passing that, in Herzfeld (1955, p. 691, equa- 
tions 18-3), the last term in this brace is cited incorrectly by omission of the 
factor f,(1—f,) which makes that term vanish for a pure gas. 

Besides giving a check on the correctness of (16), (17) and (18), equations (28) 
and (30), when boldly extrapolated to finite shock strengths, give a very useful 
approximate survey of the results which are obtained more accurately by 
numerical integrations. 


is 


10. Sample numerical calculations 

The dozen examples shown here were picked to exhibit the following. 

(a) The effect of shock strength on the concentration-versus-velocity integrals 
for a single initial gas mixture (argon—helium, equal parts by moles). The Mach 
numbers chosen were 1-34 (¢ = 1), 2-05 (¢ = 2), 4-84 (e€ = 2-8) and oo (e = 3), the 
last obviously just for fun. These results are shown in figure 2, where they are 
compared with the approximation given by equations (24) and (28). As might 
have been expected, the weak-shock approximation is quite accurate at ¢ = 1, 
and discrepancies between it and the numerical results increase continually with 
increasing ¢. The merit of the series solutions in terms of the parameter e, which 
remains finite as .W - 00, is seen from the fair approximation given by its first 
term even at very large shock strengths. 

(6) The effect of initial mixture ratio on various shock properties, for a parti- 
cular gas-pair (argon—helium) and shock strength (¢ = 2). Thus figure 3 shows 
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concentration-velocity integrals on a semi-logarithmic scale (to emphasize per- 
centage changes in concentration). Note that these are decidedly more sym- 
metrical about the mid-velocity, and thus agree better qualitatively with the 
weak-shock theory when f, is small, than when it is large. This is a coincidence 
brought about by the superposition of two effects, the first being the general 
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movement of the maximum separation towards the upstream side of the shock 
(towards the right in figure 3) in response to ordinary diffusion and the fact that 
the separating potential is Inp rather than p (as discussed in §2). The second 
effect is intuitively rather obvious, and well illustrated by figures 4 and 5. When 
the initial mixture contains just a little of the lighter gas, the separation is brought 
about by the light molecules of the shocked gas running out ahead of the main 
shock or mean-velocity transition of the entire mixture. This evidently shifts the 
maximum separation even farther upstream relative to the main shock and leads 
to the very asymmetrical f-v curves in figure 3a. When the initial mixture con- 
tains just a little heavy gas, the heavier molecules lag behind the main shock, thus 
shifting the maximum separation downstream and partially cancelling the first 
effect mentioned above. 

If the reader has been distracted from the foregoing argument by the rather 
odd appearance of the argon velocity profile in figure 5, he is forgiven forthwith. 
The behaviour of that curve, which indicates that the argon is momentarily 
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accelerated, rather than slowed down, upon entry into the shock, is the single most 
intriguing—and perhaps unbelievable—result of the present study. It implies, 
upon transformation of co-ordinates to a frame at rest in the unshocked gas, that 
upon arrival of the shock the argon first jumps backwards, then turns and runs 
fast to catch up with the shock. Thesame prediction is inherent in Dyakov’s weak- 
shock theory, which incidentally predicts a corresponding undershoot of the 
helium velocity in figure 4. It is interesting to note that the latter prediction is one 
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FicurE 4. Shock profiles with very little light gas; f, = 0-9, « = 2. 
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qualitative implication of the weak-shock theory which is refuted by the accurate 
numerical integrations. Regarding the physical plausibility of this effect, the 
author has no very fixed opinion as yet. It is certainly hard to imagine the source 
of the effect from a molecular viewpoint and no corresponding effect is found in 
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FIGURE 7. Entropy-velocity integrals; ¢ = 2. 


Carrier’s dusty-gas analysis, if the latter may be considered to give an indication 
of proper physical behaviour in a limit of very large mass ratio. A mass ratio of 
10 (argon—helium) and shock Mach number of 2:05 may simply carry the present 
analysis beyond the valid range of the Navier-Stokes approximation although 
that approximation would still be quite accurate for a Mach 2-05 shock in a pure 
monatomic gas (cf. Talbot & Sherman 1959). 














e 


n 








Shock-wave structure in binary mixtures 477 


Figures 4, 5 and 6 (the last showing velocity profiles for the intermediate case 
fy = 0-5) emphasize that a shock wave in a mixture may be quite complicated 
compared to that in a pure gas, and that it may actually extend over a distance 
which is quite large compared with the conventionally defined shock thickness 
(cf. equation (29)). They also suggest that experimentalists interested in this 
problem should take care to know exactly to what flow properties their instru- 
ments respond. 

Figure 7 shows the effect of initial mixture ratio on the entropy-velocity relation- 
ship within the shock. Actually, the quantity labelled s here is not strictly the 
entropy of the actual non-equilibrium state of the gas at a point inside the shock, 
but the entropy of the equilibrium state having the same pressure, temperature 
and concentration. This is a quantity easy to calculate—in our case from the 


S-S__, (T\Un-D y, ey 
PR = fia] (7 ) : |+a -f)\0] (7 ral (31) 


It is presumably slightly larger than the entropy of the true non-equilibrium 
state, since a bit of fluid, if isolated in the latter condition, would spontaneously 
tend to the condition whose entropy we calculate. 

Setting aside this rather fine point, we see in figure 7 just what we expect. While 
the mixture is being separated in the upstream portion of the shock, entropy rises 
at a rate smaller than in a pure gas shock of the same strength. While remixing is 
taking place farther downstream, entropy rises faster (or falls more slowly) than 
in the pure gas. 

Figure 8 shows some over-all or extreme properties of the shock transition as 
functions of f,. From 8a we see both that diffusion effects broaden a shock in 
A-He quite spectacularly, and that the weak-shock theory, equation (30), gives 
a very good idea of the dependence of this broadening upon initial mixture ratio. 
From 8b we see that the maximum ratio of diffusion velocity to mixture velocity, 
(u, —U,)/u, is smallest when diffusive broadening of the shock is largest, and that 
the highest diffusion velocities appear when there is only a trace of one gas or the 
other present. The quantity plotted is simply related to f by use of equations (5), 


(6) and (7), i.e. Uy — Us fyo-f 14+(0-1)f 


formula 





_ Jo. 5 (32) 
a fU-f)1+(O-Vfo 
The weak-shock approximation to this result is 
( M,-M, y-1 
. y+ (yy? | ” iil Sig. | 
U,— Us 7 27 + 1) (9 ) iH ) (1 — w*) 
u , Sy Se | Rs l a Vfol—fo) (-M, y-1 . ) 
i Pr Se | M, y (33) 


This equation gives a fairly true image of the numerical results, particularly for 
small f,, but is considerably worse at large f, for reasons reviewed at the beginning 
of this subsection. 

(c) The effect of mass ratios (and accompanying mixture properties) at con- 
stant f, and shock strength. Instead of making an artificial parameter study by 
changing @ alone, we deemed it more interesting to treat real mixtures of various 
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gases, to attain a wide range of 6 along with the associated changes in «, Pr, and Sc 
which nature provides. Thus calculations were made for f, = 0-5 and € = 2 in the 
mixtures #9A—36A (9 = 10/9), A-Ne (@ = 2), A-He (9 = 10) and Xe—He (6 = 32:8). 
Table 1 shows the variation of shock thickness and maximum diffusion velocity 
with @ in these mixtures, all of which have y = 5/3 and X = 4/3. 
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M,—-M, Numerical Equation Numerical F.quation | 
Gas pair 0 M, Ly Pr Se integration 49 integration 52 
Pure argon l 0 0 0-704 8-15 8°55 

10 A 36 A 10/9 0-10526 0-013 0-704 0-784 8-18 8:58 00209 0-0198 | 

A-Ne 2 0-66667 0-17 0-667 0°750 8-99 9-49 0-116 0-112 

A-—He 10 16364 037 0-440 0-330 19-84 21-15 0-273 0-284 
Xe-He 32:8 11-8817 0:434 0-250 0-200 35:8 37:8 0-276 0-301! 
TaBLE 1. Effect of @ on shock thickness and maximum diffusion velocity. fy = 0-5, MH = 2-05, y = 5/3. 

X = 4/3 for all cases; uy, = viscosity at point of maximum yp du/dx. 


We see again that the weak-shock theory is very useful for rough predictions 
at finite shock strength. Note that quite refined experimentation would be 
required to detect either the minute separation of the isotopic mixture *A—*A, 
or the broadening of the shock due to diffusion in this case. A similar result is 
found for shocks in atmospheric air, in which the importance of 6 is emphasized 
by the fact that the diffusive shock broadening due to the 0-94 % of argon present 
is about 80° as much as that due to the diffusion of the oxygen and nitrogen. 
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11. Criticism of the theory when @ is large 

As mentioned several times above, the present theory seems deficient in the 
case of large 0, when it insists that both mixture components always share the 
same temperature. This has led to some intuitively implausible results in the 
examples worked out here for argon-helium mixtures. For example, figure 4 
shows the mean-velocity change of the helium almost half-way accomplished 
before the random velocity (temperature) of this constituent has been appreciably 
influenced. One might intuitively expect these two quantities to start changing 
at about the same time. 

To get a better picture of actual events during shock passage through a mixture 
with large #, one might try at least two theoretical approaches. The first is to 
devise some supplementary ‘relaxation’ equation to account for the difficulty of 
maintaining thermal equilibrium between light and heavy constituents, and to 
couple this equation to the present set, modifying the latter somewhat to admit 
the existence of two different temperatures. An example of this approach is given 
by Jukes (1957) who discusses shock structure in a mixture of protons and 
electrons. Jukes’s calculation shows a significant difference between electron and 
ion temperatures throughout a wide zone ahead of the ‘velocity shock’ but sheds 
no direct light on the present problem, since diffusive separation of electrons and 
ions was taken to be negligible in view of the intense electric field which this would 
produce. Regarding this technique it may be noted pessimistically that the 
coupling of the ordinary sort of relaxation equation (cf. Jukes 1957) to the 
Navier-Stokes equations leads to very nasty singular-point behaviour of the 
resulting equations, which precludes straightforward numerical integration from 
one side of the shock to the other. 

The second possibly profitable approach would be an extension of the method 
of Mott-Smith (1951) to treat a gas mixture. Again, this approach has been 
applied to a hydrogen plasma (Tidman 1958), but not to a mixture in which 
diffusive separation is unrestrained by electric fields. 

Finally, of course, one might attempt an experimental determination of the 
true state of affairs with shocks in, say, A-He mixtures. In fact, some early work 
along this line, by Professor D. F. Hornig of Brown University (now of Princeton) 


gave the author his first inclination to study this problem. 
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